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 /* More checks are in triple check (grompp.c) */
1000 if (ir->coulombtype == eelSWITCH)
1002 sprintf(warn_buf, "coulombtype = %s is only for testing purposes and can lead to serious "
1003 "artifacts, advice: use coulombtype = %s",
1004 eel_names[ir->coulombtype],
1005 eel_names[eelRF_ZERO]);
1006 warning(wi, warn_buf);
1009 if (ir->epsilon_r != 1 && ir->implicit_solvent == eisGBSA)
1011 sprintf(warn_buf, "epsilon-r = %g with GB implicit solvent, will use this value for inner dielectric", ir->epsilon_r);
1012 warning_note(wi, warn_buf);
1015 if (EEL_RF(ir->coulombtype) && ir->epsilon_rf == 1 && ir->epsilon_r != 1)
1017 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);
1018 warning(wi, warn_buf);
1019 ir->epsilon_rf = ir->epsilon_r;
1020 ir->epsilon_r = 1.0;
1023 if (getenv("GALACTIC_DYNAMICS") == NULL)
1025 sprintf(err_buf, "epsilon-r must be >= 0 instead of %g\n", ir->epsilon_r);
1026 CHECK(ir->epsilon_r < 0);
1029 if (EEL_RF(ir->coulombtype))
1031 /* reaction field (at the cut-off) */
1033 if (ir->coulombtype == eelRF_ZERO)
1035 sprintf(warn_buf, "With coulombtype = %s, epsilon-rf must be 0, assuming you meant epsilon_rf=0",
1036 eel_names[ir->coulombtype]);
1037 CHECK(ir->epsilon_rf != 0);
1038 ir->epsilon_rf = 0.0;
1041 sprintf(err_buf, "epsilon-rf must be >= epsilon-r");
1042 CHECK((ir->epsilon_rf < ir->epsilon_r && ir->epsilon_rf != 0) ||
1043 (ir->epsilon_r == 0));
1044 if (ir->epsilon_rf == ir->epsilon_r)
1046 sprintf(warn_buf, "Using epsilon-rf = epsilon-r with %s does not make sense",
1047 eel_names[ir->coulombtype]);
1048 warning(wi, warn_buf);
1051 /* Allow rlist>rcoulomb for tabulated long range stuff. This just
1052 * means the interaction is zero outside rcoulomb, but it helps to
1053 * provide accurate energy conservation.
1055 if (EEL_MIGHT_BE_ZERO_AT_CUTOFF(ir->coulombtype))
1057 if (EEL_SWITCHED(ir->coulombtype))
1060 "With coulombtype = %s rcoulomb_switch must be < rcoulomb. Or, better: Use the potential modifier options!",
1061 eel_names[ir->coulombtype]);
1062 CHECK(ir->rcoulomb_switch >= ir->rcoulomb);
1065 else if (ir->coulombtype == eelCUT || EEL_RF(ir->coulombtype))
1067 if (ir->cutoff_scheme == ecutsGROUP && ir->coulomb_modifier == eintmodNONE)
1069 sprintf(err_buf, "With coulombtype = %s, rcoulomb should be >= rlist unless you use a potential modifier",
1070 eel_names[ir->coulombtype]);
1071 CHECK(ir->rlist > ir->rcoulomb);
1075 if (ir->coulombtype == eelSWITCH || ir->coulombtype == eelSHIFT ||
1076 ir->vdwtype == evdwSWITCH || ir->vdwtype == evdwSHIFT)
1079 "The switch/shift interaction settings are just for compatibility; you will get better "
1080 "performance from applying potential modifiers to your interactions!\n");
1081 warning_note(wi, warn_buf);
1084 if (ir->coulombtype == eelPMESWITCH)
1086 if (ir->rcoulomb_switch/ir->rcoulomb < 0.9499)
1088 sprintf(warn_buf, "The switching range for %s should be 5%% or less, energy conservation will be good anyhow, since ewald_rtol = %g",
1089 eel_names[ir->coulombtype],
1091 warning(wi, warn_buf);
1095 if (EEL_FULL(ir->coulombtype))
1097 if (ir->coulombtype == eelPMESWITCH || ir->coulombtype == eelPMEUSER ||
1098 ir->coulombtype == eelPMEUSERSWITCH)
1100 sprintf(err_buf, "With coulombtype = %s, rcoulomb must be <= rlist",
1101 eel_names[ir->coulombtype]);
1102 CHECK(ir->rcoulomb > ir->rlist);
1104 else if (ir->cutoff_scheme == ecutsGROUP && ir->coulomb_modifier == eintmodNONE)
1106 if (ir->coulombtype == eelPME || ir->coulombtype == eelP3M_AD)
1109 "With coulombtype = %s (without modifier), rcoulomb must be equal to rlist,\n"
1110 "or rlistlong if nstcalclr=1. For optimal energy conservation,consider using\n"
1111 "a potential modifier.", eel_names[ir->coulombtype]);
1112 if (ir->nstcalclr == 1)
1114 CHECK(ir->rcoulomb != ir->rlist && ir->rcoulomb != ir->rlistlong);
1118 CHECK(ir->rcoulomb != ir->rlist);
1124 if (EEL_PME(ir->coulombtype))
1126 if (ir->pme_order < 3)
1128 warning_error(wi, "pme-order can not be smaller than 3");
1132 if (ir->nwall == 2 && EEL_FULL(ir->coulombtype))
1134 if (ir->ewald_geometry == eewg3D)
1136 sprintf(warn_buf, "With pbc=%s you should use ewald-geometry=%s",
1137 epbc_names[ir->ePBC], eewg_names[eewg3DC]);
1138 warning(wi, warn_buf);
1140 /* This check avoids extra pbc coding for exclusion corrections */
1141 sprintf(err_buf, "wall-ewald-zfac should be >= 2");
1142 CHECK(ir->wall_ewald_zfac < 2);
1145 if (EVDW_SWITCHED(ir->vdwtype))
1147 sprintf(err_buf, "With vdwtype = %s rvdw-switch must be < rvdw. Or, better - use a potential modifier.",
1148 evdw_names[ir->vdwtype]);
1149 CHECK(ir->rvdw_switch >= ir->rvdw);
1151 else if (ir->vdwtype == evdwCUT)
1153 if (ir->cutoff_scheme == ecutsGROUP && ir->vdw_modifier == eintmodNONE)
1155 sprintf(err_buf, "With vdwtype = %s, rvdw must be >= rlist unless you use a potential modifier", evdw_names[ir->vdwtype]);
1156 CHECK(ir->rlist > ir->rvdw);
1159 if (ir->cutoff_scheme == ecutsGROUP)
1161 if (((ir->coulomb_modifier != eintmodNONE && ir->rcoulomb == ir->rlist) ||
1162 (ir->vdw_modifier != eintmodNONE && ir->rvdw == ir->rlist)) &&
1165 warning_note(wi, "With exact cut-offs, rlist should be "
1166 "larger than rcoulomb and rvdw, so that there "
1167 "is a buffer region for particle motion "
1168 "between neighborsearch steps");
1171 if (EEL_IS_ZERO_AT_CUTOFF(ir->coulombtype)
1172 && (ir->rlistlong <= ir->rcoulomb))
1174 sprintf(warn_buf, "For energy conservation with switch/shift potentials, %s should be 0.1 to 0.3 nm larger than rcoulomb.",
1175 IR_TWINRANGE(*ir) ? "rlistlong" : "rlist");
1176 warning_note(wi, warn_buf);
1178 if (EVDW_SWITCHED(ir->vdwtype) && (ir->rlistlong <= ir->rvdw))
1180 sprintf(warn_buf, "For energy conservation with switch/shift potentials, %s should be 0.1 to 0.3 nm larger than rvdw.",
1181 IR_TWINRANGE(*ir) ? "rlistlong" : "rlist");
1182 warning_note(wi, warn_buf);
1186 if (ir->vdwtype == evdwUSER && ir->eDispCorr != edispcNO)
1188 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.");
1191 if (ir->nstlist == -1)
1193 sprintf(err_buf, "With nstlist=-1 rvdw and rcoulomb should be smaller than rlist to account for diffusion and possibly charge-group radii");
1194 CHECK(ir->rvdw >= ir->rlist || ir->rcoulomb >= ir->rlist);
1196 sprintf(err_buf, "nstlist can not be smaller than -1");
1197 CHECK(ir->nstlist < -1);
1199 if (ir->eI == eiLBFGS && (ir->coulombtype == eelCUT || ir->vdwtype == evdwCUT)
1202 warning(wi, "For efficient BFGS minimization, use switch/shift/pme instead of cut-off.");
1205 if (ir->eI == eiLBFGS && ir->nbfgscorr <= 0)
1207 warning(wi, "Using L-BFGS with nbfgscorr<=0 just gets you steepest descent.");
1210 /* ENERGY CONSERVATION */
1211 if (ir_NVE(ir) && ir->cutoff_scheme == ecutsGROUP)
1213 if (!EVDW_MIGHT_BE_ZERO_AT_CUTOFF(ir->vdwtype) && ir->rvdw > 0 && ir->vdw_modifier == eintmodNONE)
1215 sprintf(warn_buf, "You are using a cut-off for VdW interactions with NVE, for good energy conservation use vdwtype = %s (possibly with DispCorr)",
1216 evdw_names[evdwSHIFT]);
1217 warning_note(wi, warn_buf);
1219 if (!EEL_MIGHT_BE_ZERO_AT_CUTOFF(ir->coulombtype) && ir->rcoulomb > 0 && ir->coulomb_modifier == eintmodNONE)
1221 sprintf(warn_buf, "You are using a cut-off for electrostatics with NVE, for good energy conservation use coulombtype = %s or %s",
1222 eel_names[eelPMESWITCH], eel_names[eelRF_ZERO]);
1223 warning_note(wi, warn_buf);
1227 /* IMPLICIT SOLVENT */
1228 if (ir->coulombtype == eelGB_NOTUSED)
1230 ir->coulombtype = eelCUT;
1231 ir->implicit_solvent = eisGBSA;
1232 fprintf(stderr, "Note: Old option for generalized born electrostatics given:\n"
1233 "Changing coulombtype from \"generalized-born\" to \"cut-off\" and instead\n"
1234 "setting implicit-solvent value to \"GBSA\" in input section.\n");
1237 if (ir->sa_algorithm == esaSTILL)
1239 sprintf(err_buf, "Still SA algorithm not available yet, use %s or %s instead\n", esa_names[esaAPPROX], esa_names[esaNO]);
1240 CHECK(ir->sa_algorithm == esaSTILL);
1243 if (ir->implicit_solvent == eisGBSA)
1245 sprintf(err_buf, "With GBSA implicit solvent, rgbradii must be equal to rlist.");
1246 CHECK(ir->rgbradii != ir->rlist);
1248 if (ir->coulombtype != eelCUT)
1250 sprintf(err_buf, "With GBSA, coulombtype must be equal to %s\n", eel_names[eelCUT]);
1251 CHECK(ir->coulombtype != eelCUT);
1253 if (ir->vdwtype != evdwCUT)
1255 sprintf(err_buf, "With GBSA, vdw-type must be equal to %s\n", evdw_names[evdwCUT]);
1256 CHECK(ir->vdwtype != evdwCUT);
1258 if (ir->nstgbradii < 1)
1260 sprintf(warn_buf, "Using GBSA with nstgbradii<1, setting nstgbradii=1");
1261 warning_note(wi, warn_buf);
1264 if (ir->sa_algorithm == esaNO)
1266 sprintf(warn_buf, "No SA (non-polar) calculation requested together with GB. Are you sure this is what you want?\n");
1267 warning_note(wi, warn_buf);
1269 if (ir->sa_surface_tension < 0 && ir->sa_algorithm != esaNO)
1271 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");
1272 warning_note(wi, warn_buf);
1274 if (ir->gb_algorithm == egbSTILL)
1276 ir->sa_surface_tension = 0.0049 * CAL2JOULE * 100;
1280 ir->sa_surface_tension = 0.0054 * CAL2JOULE * 100;
1283 if (ir->sa_surface_tension == 0 && ir->sa_algorithm != esaNO)
1285 sprintf(err_buf, "Surface tension set to 0 while SA-calculation requested\n");
1286 CHECK(ir->sa_surface_tension == 0 && ir->sa_algorithm != esaNO);
1293 if (ir->cutoff_scheme != ecutsGROUP)
1295 warning_error(wi, "AdresS simulation supports only cutoff-scheme=group");
1299 warning_error(wi, "AdresS simulation supports only stochastic dynamics");
1301 if (ir->epc != epcNO)
1303 warning_error(wi, "AdresS simulation does not support pressure coupling");
1305 if (EEL_FULL(ir->coulombtype))
1307 warning_error(wi, "AdresS simulation does not support long-range electrostatics");
1312 /* count the number of text elemets separated by whitespace in a string.
1313 str = the input string
1314 maxptr = the maximum number of allowed elements
1315 ptr = the output array of pointers to the first character of each element
1316 returns: the number of elements. */
1317 int str_nelem(const char *str, int maxptr, char *ptr[])
1322 copy0 = strdup(str);
1325 while (*copy != '\0')
1329 gmx_fatal(FARGS, "Too many groups on line: '%s' (max is %d)",
1337 while ((*copy != '\0') && !isspace(*copy))
1356 /* interpret a number of doubles from a string and put them in an array,
1357 after allocating space for them.
1358 str = the input string
1359 n = the (pre-allocated) number of doubles read
1360 r = the output array of doubles. */
1361 static void parse_n_real(char *str, int *n, real **r)
1366 *n = str_nelem(str, MAXPTR, ptr);
1369 for (i = 0; i < *n; i++)
1371 (*r)[i] = strtod(ptr[i], NULL);
1375 static void do_fep_params(t_inputrec *ir, char fep_lambda[][STRLEN], char weights[STRLEN])
1378 int i, j, max_n_lambda, nweights, nfep[efptNR];
1379 t_lambda *fep = ir->fepvals;
1380 t_expanded *expand = ir->expandedvals;
1381 real **count_fep_lambdas;
1382 gmx_bool bOneLambda = TRUE;
1384 snew(count_fep_lambdas, efptNR);
1386 /* FEP input processing */
1387 /* first, identify the number of lambda values for each type.
1388 All that are nonzero must have the same number */
1390 for (i = 0; i < efptNR; i++)
1392 parse_n_real(fep_lambda[i], &(nfep[i]), &(count_fep_lambdas[i]));
1395 /* now, determine the number of components. All must be either zero, or equal. */
1398 for (i = 0; i < efptNR; i++)
1400 if (nfep[i] > max_n_lambda)
1402 max_n_lambda = nfep[i]; /* here's a nonzero one. All of them
1403 must have the same number if its not zero.*/
1408 for (i = 0; i < efptNR; i++)
1412 ir->fepvals->separate_dvdl[i] = FALSE;
1414 else if (nfep[i] == max_n_lambda)
1416 if (i != efptTEMPERATURE) /* we treat this differently -- not really a reason to compute the derivative with
1417 respect to the temperature currently */
1419 ir->fepvals->separate_dvdl[i] = TRUE;
1424 gmx_fatal(FARGS, "Number of lambdas (%d) for FEP type %s not equal to number of other types (%d)",
1425 nfep[i], efpt_names[i], max_n_lambda);
1428 /* we don't print out dhdl if the temperature is changing, since we can't correctly define dhdl in this case */
1429 ir->fepvals->separate_dvdl[efptTEMPERATURE] = FALSE;
1431 /* the number of lambdas is the number we've read in, which is either zero
1432 or the same for all */
1433 fep->n_lambda = max_n_lambda;
1435 /* allocate space for the array of lambda values */
1436 snew(fep->all_lambda, efptNR);
1437 /* if init_lambda is defined, we need to set lambda */
1438 if ((fep->init_lambda > 0) && (fep->n_lambda == 0))
1440 ir->fepvals->separate_dvdl[efptFEP] = TRUE;
1442 /* otherwise allocate the space for all of the lambdas, and transfer the data */
1443 for (i = 0; i < efptNR; i++)
1445 snew(fep->all_lambda[i], fep->n_lambda);
1446 if (nfep[i] > 0) /* if it's zero, then the count_fep_lambda arrays
1449 for (j = 0; j < fep->n_lambda; j++)
1451 fep->all_lambda[i][j] = (double)count_fep_lambdas[i][j];
1453 sfree(count_fep_lambdas[i]);
1456 sfree(count_fep_lambdas);
1458 /* "fep-vals" is either zero or the full number. If zero, we'll need to define fep-lambdas for internal
1459 bookkeeping -- for now, init_lambda */
1461 if ((nfep[efptFEP] == 0) && (fep->init_lambda >= 0))
1463 for (i = 0; i < fep->n_lambda; i++)
1465 fep->all_lambda[efptFEP][i] = fep->init_lambda;
1469 /* check to see if only a single component lambda is defined, and soft core is defined.
1470 In this case, turn on coulomb soft core */
1472 if (max_n_lambda == 0)
1478 for (i = 0; i < efptNR; i++)
1480 if ((nfep[i] != 0) && (i != efptFEP))
1486 if ((bOneLambda) && (fep->sc_alpha > 0))
1488 fep->bScCoul = TRUE;
1491 /* Fill in the others with the efptFEP if they are not explicitly
1492 specified (i.e. nfep[i] == 0). This means if fep is not defined,
1493 they are all zero. */
1495 for (i = 0; i < efptNR; i++)
1497 if ((nfep[i] == 0) && (i != efptFEP))
1499 for (j = 0; j < fep->n_lambda; j++)
1501 fep->all_lambda[i][j] = fep->all_lambda[efptFEP][j];
1507 /* make it easier if sc_r_power = 48 by increasing it to the 4th power, to be in the right scale. */
1508 if (fep->sc_r_power == 48)
1510 if (fep->sc_alpha > 0.1)
1512 gmx_fatal(FARGS, "sc_alpha (%f) for sc_r_power = 48 should usually be between 0.001 and 0.004", fep->sc_alpha);
1516 expand = ir->expandedvals;
1517 /* now read in the weights */
1518 parse_n_real(weights, &nweights, &(expand->init_lambda_weights));
1521 snew(expand->init_lambda_weights, fep->n_lambda); /* initialize to zero */
1523 else if (nweights != fep->n_lambda)
1525 gmx_fatal(FARGS, "Number of weights (%d) is not equal to number of lambda values (%d)",
1526 nweights, fep->n_lambda);
1528 if ((expand->nstexpanded < 0) && (ir->efep != efepNO))
1530 expand->nstexpanded = fep->nstdhdl;
1531 /* if you don't specify nstexpanded when doing expanded ensemble free energy calcs, it is set to nstdhdl */
1533 if ((expand->nstexpanded < 0) && ir->bSimTemp)
1535 expand->nstexpanded = 2*(int)(ir->opts.tau_t[0]/ir->delta_t);
1536 /* if you don't specify nstexpanded when doing expanded ensemble simulated tempering, it is set to
1537 2*tau_t just to be careful so it's not to frequent */
1542 static void do_simtemp_params(t_inputrec *ir)
1545 snew(ir->simtempvals->temperatures, ir->fepvals->n_lambda);
1546 GetSimTemps(ir->fepvals->n_lambda, ir->simtempvals, ir->fepvals->all_lambda[efptTEMPERATURE]);
1551 static void do_wall_params(t_inputrec *ir,
1552 char *wall_atomtype, char *wall_density,
1556 char *names[MAXPTR];
1559 opts->wall_atomtype[0] = NULL;
1560 opts->wall_atomtype[1] = NULL;
1562 ir->wall_atomtype[0] = -1;
1563 ir->wall_atomtype[1] = -1;
1564 ir->wall_density[0] = 0;
1565 ir->wall_density[1] = 0;
1569 nstr = str_nelem(wall_atomtype, MAXPTR, names);
1570 if (nstr != ir->nwall)
1572 gmx_fatal(FARGS, "Expected %d elements for wall_atomtype, found %d",
1575 for (i = 0; i < ir->nwall; i++)
1577 opts->wall_atomtype[i] = strdup(names[i]);
1580 if (ir->wall_type == ewt93 || ir->wall_type == ewt104)
1582 nstr = str_nelem(wall_density, MAXPTR, names);
1583 if (nstr != ir->nwall)
1585 gmx_fatal(FARGS, "Expected %d elements for wall-density, found %d", ir->nwall, nstr);
1587 for (i = 0; i < ir->nwall; i++)
1589 sscanf(names[i], "%lf", &dbl);
1592 gmx_fatal(FARGS, "wall-density[%d] = %f\n", i, dbl);
1594 ir->wall_density[i] = dbl;
1600 static void add_wall_energrps(gmx_groups_t *groups, int nwall, t_symtab *symtab)
1608 srenew(groups->grpname, groups->ngrpname+nwall);
1609 grps = &(groups->grps[egcENER]);
1610 srenew(grps->nm_ind, grps->nr+nwall);
1611 for (i = 0; i < nwall; i++)
1613 sprintf(str, "wall%d", i);
1614 groups->grpname[groups->ngrpname] = put_symtab(symtab, str);
1615 grps->nm_ind[grps->nr++] = groups->ngrpname++;
1620 void read_expandedparams(int *ninp_p, t_inpfile **inp_p,
1621 t_expanded *expand, warninp_t wi)
1623 int ninp, nerror = 0;
1629 /* read expanded ensemble parameters */
1630 CCTYPE ("expanded ensemble variables");
1631 ITYPE ("nstexpanded", expand->nstexpanded, -1);
1632 EETYPE("lmc-stats", expand->elamstats, elamstats_names);
1633 EETYPE("lmc-move", expand->elmcmove, elmcmove_names);
1634 EETYPE("lmc-weights-equil", expand->elmceq, elmceq_names);
1635 ITYPE ("weight-equil-number-all-lambda", expand->equil_n_at_lam, -1);
1636 ITYPE ("weight-equil-number-samples", expand->equil_samples, -1);
1637 ITYPE ("weight-equil-number-steps", expand->equil_steps, -1);
1638 RTYPE ("weight-equil-wl-delta", expand->equil_wl_delta, -1);
1639 RTYPE ("weight-equil-count-ratio", expand->equil_ratio, -1);
1640 CCTYPE("Seed for Monte Carlo in lambda space");
1641 ITYPE ("lmc-seed", expand->lmc_seed, -1);
1642 RTYPE ("mc-temperature", expand->mc_temp, -1);
1643 ITYPE ("lmc-repeats", expand->lmc_repeats, 1);
1644 ITYPE ("lmc-gibbsdelta", expand->gibbsdeltalam, -1);
1645 ITYPE ("lmc-forced-nstart", expand->lmc_forced_nstart, 0);
1646 EETYPE("symmetrized-transition-matrix", expand->bSymmetrizedTMatrix, yesno_names);
1647 ITYPE("nst-transition-matrix", expand->nstTij, -1);
1648 ITYPE ("mininum-var-min", expand->minvarmin, 100); /*default is reasonable */
1649 ITYPE ("weight-c-range", expand->c_range, 0); /* default is just C=0 */
1650 RTYPE ("wl-scale", expand->wl_scale, 0.8);
1651 RTYPE ("wl-ratio", expand->wl_ratio, 0.8);
1652 RTYPE ("init-wl-delta", expand->init_wl_delta, 1.0);
1653 EETYPE("wl-oneovert", expand->bWLoneovert, yesno_names);
1661 void get_ir(const char *mdparin, const char *mdparout,
1662 t_inputrec *ir, t_gromppopts *opts,
1666 double dumdub[2][6];
1670 char warn_buf[STRLEN];
1671 t_lambda *fep = ir->fepvals;
1672 t_expanded *expand = ir->expandedvals;
1674 inp = read_inpfile(mdparin, &ninp, NULL, wi);
1676 snew(dumstr[0], STRLEN);
1677 snew(dumstr[1], STRLEN);
1679 /* remove the following deprecated commands */
1682 REM_TYPE("domain-decomposition");
1683 REM_TYPE("andersen-seed");
1685 REM_TYPE("dihre-fc");
1686 REM_TYPE("dihre-tau");
1687 REM_TYPE("nstdihreout");
1688 REM_TYPE("nstcheckpoint");
1690 /* replace the following commands with the clearer new versions*/
1691 REPL_TYPE("unconstrained-start", "continuation");
1692 REPL_TYPE("foreign-lambda", "fep-lambdas");
1694 CCTYPE ("VARIOUS PREPROCESSING OPTIONS");
1695 CTYPE ("Preprocessor information: use cpp syntax.");
1696 CTYPE ("e.g.: -I/home/joe/doe -I/home/mary/roe");
1697 STYPE ("include", opts->include, NULL);
1698 CTYPE ("e.g.: -DPOSRES -DFLEXIBLE (note these variable names are case sensitive)");
1699 STYPE ("define", opts->define, NULL);
1701 CCTYPE ("RUN CONTROL PARAMETERS");
1702 EETYPE("integrator", ir->eI, ei_names);
1703 CTYPE ("Start time and timestep in ps");
1704 RTYPE ("tinit", ir->init_t, 0.0);
1705 RTYPE ("dt", ir->delta_t, 0.001);
1706 STEPTYPE ("nsteps", ir->nsteps, 0);
1707 CTYPE ("For exact run continuation or redoing part of a run");
1708 STEPTYPE ("init-step", ir->init_step, 0);
1709 CTYPE ("Part index is updated automatically on checkpointing (keeps files separate)");
1710 ITYPE ("simulation-part", ir->simulation_part, 1);
1711 CTYPE ("mode for center of mass motion removal");
1712 EETYPE("comm-mode", ir->comm_mode, ecm_names);
1713 CTYPE ("number of steps for center of mass motion removal");
1714 ITYPE ("nstcomm", ir->nstcomm, 100);
1715 CTYPE ("group(s) for center of mass motion removal");
1716 STYPE ("comm-grps", vcm, NULL);
1718 CCTYPE ("LANGEVIN DYNAMICS OPTIONS");
1719 CTYPE ("Friction coefficient (amu/ps) and random seed");
1720 RTYPE ("bd-fric", ir->bd_fric, 0.0);
1721 ITYPE ("ld-seed", ir->ld_seed, 1993);
1724 CCTYPE ("ENERGY MINIMIZATION OPTIONS");
1725 CTYPE ("Force tolerance and initial step-size");
1726 RTYPE ("emtol", ir->em_tol, 10.0);
1727 RTYPE ("emstep", ir->em_stepsize, 0.01);
1728 CTYPE ("Max number of iterations in relax-shells");
1729 ITYPE ("niter", ir->niter, 20);
1730 CTYPE ("Step size (ps^2) for minimization of flexible constraints");
1731 RTYPE ("fcstep", ir->fc_stepsize, 0);
1732 CTYPE ("Frequency of steepest descents steps when doing CG");
1733 ITYPE ("nstcgsteep", ir->nstcgsteep, 1000);
1734 ITYPE ("nbfgscorr", ir->nbfgscorr, 10);
1736 CCTYPE ("TEST PARTICLE INSERTION OPTIONS");
1737 RTYPE ("rtpi", ir->rtpi, 0.05);
1739 /* Output options */
1740 CCTYPE ("OUTPUT CONTROL OPTIONS");
1741 CTYPE ("Output frequency for coords (x), velocities (v) and forces (f)");
1742 ITYPE ("nstxout", ir->nstxout, 0);
1743 ITYPE ("nstvout", ir->nstvout, 0);
1744 ITYPE ("nstfout", ir->nstfout, 0);
1745 ir->nstcheckpoint = 1000;
1746 CTYPE ("Output frequency for energies to log file and energy file");
1747 ITYPE ("nstlog", ir->nstlog, 1000);
1748 ITYPE ("nstcalcenergy", ir->nstcalcenergy, 100);
1749 ITYPE ("nstenergy", ir->nstenergy, 1000);
1750 CTYPE ("Output frequency and precision for .xtc file");
1751 ITYPE ("nstxtcout", ir->nstxtcout, 0);
1752 RTYPE ("xtc-precision", ir->xtcprec, 1000.0);
1753 CTYPE ("This selects the subset of atoms for the .xtc file. You can");
1754 CTYPE ("select multiple groups. By default all atoms will be written.");
1755 STYPE ("xtc-grps", xtc_grps, NULL);
1756 CTYPE ("Selection of energy groups");
1757 STYPE ("energygrps", energy, NULL);
1759 /* Neighbor searching */
1760 CCTYPE ("NEIGHBORSEARCHING PARAMETERS");
1761 CTYPE ("cut-off scheme (group: using charge groups, Verlet: particle based cut-offs)");
1762 EETYPE("cutoff-scheme", ir->cutoff_scheme, ecutscheme_names);
1763 CTYPE ("nblist update frequency");
1764 ITYPE ("nstlist", ir->nstlist, 10);
1765 CTYPE ("ns algorithm (simple or grid)");
1766 EETYPE("ns-type", ir->ns_type, ens_names);
1767 /* set ndelta to the optimal value of 2 */
1769 CTYPE ("Periodic boundary conditions: xyz, no, xy");
1770 EETYPE("pbc", ir->ePBC, epbc_names);
1771 EETYPE("periodic-molecules", ir->bPeriodicMols, yesno_names);
1772 CTYPE ("Allowed energy drift due to the Verlet buffer in kJ/mol/ps per atom,");
1773 CTYPE ("a value of -1 means: use rlist");
1774 RTYPE("verlet-buffer-drift", ir->verletbuf_drift, 0.005);
1775 CTYPE ("nblist cut-off");
1776 RTYPE ("rlist", ir->rlist, 1.0);
1777 CTYPE ("long-range cut-off for switched potentials");
1778 RTYPE ("rlistlong", ir->rlistlong, -1);
1779 ITYPE ("nstcalclr", ir->nstcalclr, -1);
1781 /* Electrostatics */
1782 CCTYPE ("OPTIONS FOR ELECTROSTATICS AND VDW");
1783 CTYPE ("Method for doing electrostatics");
1784 EETYPE("coulombtype", ir->coulombtype, eel_names);
1785 EETYPE("coulomb-modifier", ir->coulomb_modifier, eintmod_names);
1786 CTYPE ("cut-off lengths");
1787 RTYPE ("rcoulomb-switch", ir->rcoulomb_switch, 0.0);
1788 RTYPE ("rcoulomb", ir->rcoulomb, 1.0);
1789 CTYPE ("Relative dielectric constant for the medium and the reaction field");
1790 RTYPE ("epsilon-r", ir->epsilon_r, 1.0);
1791 RTYPE ("epsilon-rf", ir->epsilon_rf, 0.0);
1792 CTYPE ("Method for doing Van der Waals");
1793 EETYPE("vdw-type", ir->vdwtype, evdw_names);
1794 EETYPE("vdw-modifier", ir->vdw_modifier, eintmod_names);
1795 CTYPE ("cut-off lengths");
1796 RTYPE ("rvdw-switch", ir->rvdw_switch, 0.0);
1797 RTYPE ("rvdw", ir->rvdw, 1.0);
1798 CTYPE ("Apply long range dispersion corrections for Energy and Pressure");
1799 EETYPE("DispCorr", ir->eDispCorr, edispc_names);
1800 CTYPE ("Extension of the potential lookup tables beyond the cut-off");
1801 RTYPE ("table-extension", ir->tabext, 1.0);
1802 CTYPE ("Separate tables between energy group pairs");
1803 STYPE ("energygrp-table", egptable, NULL);
1804 CTYPE ("Spacing for the PME/PPPM FFT grid");
1805 RTYPE ("fourierspacing", ir->fourier_spacing, 0.12);
1806 CTYPE ("FFT grid size, when a value is 0 fourierspacing will be used");
1807 ITYPE ("fourier-nx", ir->nkx, 0);
1808 ITYPE ("fourier-ny", ir->nky, 0);
1809 ITYPE ("fourier-nz", ir->nkz, 0);
1810 CTYPE ("EWALD/PME/PPPM parameters");
1811 ITYPE ("pme-order", ir->pme_order, 4);
1812 RTYPE ("ewald-rtol", ir->ewald_rtol, 0.00001);
1813 EETYPE("ewald-geometry", ir->ewald_geometry, eewg_names);
1814 RTYPE ("epsilon-surface", ir->epsilon_surface, 0.0);
1815 EETYPE("optimize-fft", ir->bOptFFT, yesno_names);
1817 CCTYPE("IMPLICIT SOLVENT ALGORITHM");
1818 EETYPE("implicit-solvent", ir->implicit_solvent, eis_names);
1820 CCTYPE ("GENERALIZED BORN ELECTROSTATICS");
1821 CTYPE ("Algorithm for calculating Born radii");
1822 EETYPE("gb-algorithm", ir->gb_algorithm, egb_names);
1823 CTYPE ("Frequency of calculating the Born radii inside rlist");
1824 ITYPE ("nstgbradii", ir->nstgbradii, 1);
1825 CTYPE ("Cutoff for Born radii calculation; the contribution from atoms");
1826 CTYPE ("between rlist and rgbradii is updated every nstlist steps");
1827 RTYPE ("rgbradii", ir->rgbradii, 1.0);
1828 CTYPE ("Dielectric coefficient of the implicit solvent");
1829 RTYPE ("gb-epsilon-solvent", ir->gb_epsilon_solvent, 80.0);
1830 CTYPE ("Salt concentration in M for Generalized Born models");
1831 RTYPE ("gb-saltconc", ir->gb_saltconc, 0.0);
1832 CTYPE ("Scaling factors used in the OBC GB model. Default values are OBC(II)");
1833 RTYPE ("gb-obc-alpha", ir->gb_obc_alpha, 1.0);
1834 RTYPE ("gb-obc-beta", ir->gb_obc_beta, 0.8);
1835 RTYPE ("gb-obc-gamma", ir->gb_obc_gamma, 4.85);
1836 RTYPE ("gb-dielectric-offset", ir->gb_dielectric_offset, 0.009);
1837 EETYPE("sa-algorithm", ir->sa_algorithm, esa_names);
1838 CTYPE ("Surface tension (kJ/mol/nm^2) for the SA (nonpolar surface) part of GBSA");
1839 CTYPE ("The value -1 will set default value for Still/HCT/OBC GB-models.");
1840 RTYPE ("sa-surface-tension", ir->sa_surface_tension, -1);
1842 /* Coupling stuff */
1843 CCTYPE ("OPTIONS FOR WEAK COUPLING ALGORITHMS");
1844 CTYPE ("Temperature coupling");
1845 EETYPE("tcoupl", ir->etc, etcoupl_names);
1846 ITYPE ("nsttcouple", ir->nsttcouple, -1);
1847 ITYPE("nh-chain-length", ir->opts.nhchainlength, NHCHAINLENGTH);
1848 EETYPE("print-nose-hoover-chain-variables", ir->bPrintNHChains, yesno_names);
1849 CTYPE ("Groups to couple separately");
1850 STYPE ("tc-grps", tcgrps, NULL);
1851 CTYPE ("Time constant (ps) and reference temperature (K)");
1852 STYPE ("tau-t", tau_t, NULL);
1853 STYPE ("ref-t", ref_t, NULL);
1854 CTYPE ("pressure coupling");
1855 EETYPE("pcoupl", ir->epc, epcoupl_names);
1856 EETYPE("pcoupltype", ir->epct, epcoupltype_names);
1857 ITYPE ("nstpcouple", ir->nstpcouple, -1);
1858 CTYPE ("Time constant (ps), compressibility (1/bar) and reference P (bar)");
1859 RTYPE ("tau-p", ir->tau_p, 1.0);
1860 STYPE ("compressibility", dumstr[0], NULL);
1861 STYPE ("ref-p", dumstr[1], NULL);
1862 CTYPE ("Scaling of reference coordinates, No, All or COM");
1863 EETYPE ("refcoord-scaling", ir->refcoord_scaling, erefscaling_names);
1866 CCTYPE ("OPTIONS FOR QMMM calculations");
1867 EETYPE("QMMM", ir->bQMMM, yesno_names);
1868 CTYPE ("Groups treated Quantum Mechanically");
1869 STYPE ("QMMM-grps", QMMM, NULL);
1870 CTYPE ("QM method");
1871 STYPE("QMmethod", QMmethod, NULL);
1872 CTYPE ("QMMM scheme");
1873 EETYPE("QMMMscheme", ir->QMMMscheme, eQMMMscheme_names);
1874 CTYPE ("QM basisset");
1875 STYPE("QMbasis", QMbasis, NULL);
1876 CTYPE ("QM charge");
1877 STYPE ("QMcharge", QMcharge, NULL);
1878 CTYPE ("QM multiplicity");
1879 STYPE ("QMmult", QMmult, NULL);
1880 CTYPE ("Surface Hopping");
1881 STYPE ("SH", bSH, NULL);
1882 CTYPE ("CAS space options");
1883 STYPE ("CASorbitals", CASorbitals, NULL);
1884 STYPE ("CASelectrons", CASelectrons, NULL);
1885 STYPE ("SAon", SAon, NULL);
1886 STYPE ("SAoff", SAoff, NULL);
1887 STYPE ("SAsteps", SAsteps, NULL);
1888 CTYPE ("Scale factor for MM charges");
1889 RTYPE ("MMChargeScaleFactor", ir->scalefactor, 1.0);
1890 CTYPE ("Optimization of QM subsystem");
1891 STYPE ("bOPT", bOPT, NULL);
1892 STYPE ("bTS", bTS, NULL);
1894 /* Simulated annealing */
1895 CCTYPE("SIMULATED ANNEALING");
1896 CTYPE ("Type of annealing for each temperature group (no/single/periodic)");
1897 STYPE ("annealing", anneal, NULL);
1898 CTYPE ("Number of time points to use for specifying annealing in each group");
1899 STYPE ("annealing-npoints", anneal_npoints, NULL);
1900 CTYPE ("List of times at the annealing points for each group");
1901 STYPE ("annealing-time", anneal_time, NULL);
1902 CTYPE ("Temp. at each annealing point, for each group.");
1903 STYPE ("annealing-temp", anneal_temp, NULL);
1906 CCTYPE ("GENERATE VELOCITIES FOR STARTUP RUN");
1907 EETYPE("gen-vel", opts->bGenVel, yesno_names);
1908 RTYPE ("gen-temp", opts->tempi, 300.0);
1909 ITYPE ("gen-seed", opts->seed, 173529);
1912 CCTYPE ("OPTIONS FOR BONDS");
1913 EETYPE("constraints", opts->nshake, constraints);
1914 CTYPE ("Type of constraint algorithm");
1915 EETYPE("constraint-algorithm", ir->eConstrAlg, econstr_names);
1916 CTYPE ("Do not constrain the start configuration");
1917 EETYPE("continuation", ir->bContinuation, yesno_names);
1918 CTYPE ("Use successive overrelaxation to reduce the number of shake iterations");
1919 EETYPE("Shake-SOR", ir->bShakeSOR, yesno_names);
1920 CTYPE ("Relative tolerance of shake");
1921 RTYPE ("shake-tol", ir->shake_tol, 0.0001);
1922 CTYPE ("Highest order in the expansion of the constraint coupling matrix");
1923 ITYPE ("lincs-order", ir->nProjOrder, 4);
1924 CTYPE ("Number of iterations in the final step of LINCS. 1 is fine for");
1925 CTYPE ("normal simulations, but use 2 to conserve energy in NVE runs.");
1926 CTYPE ("For energy minimization with constraints it should be 4 to 8.");
1927 ITYPE ("lincs-iter", ir->nLincsIter, 1);
1928 CTYPE ("Lincs will write a warning to the stderr if in one step a bond");
1929 CTYPE ("rotates over more degrees than");
1930 RTYPE ("lincs-warnangle", ir->LincsWarnAngle, 30.0);
1931 CTYPE ("Convert harmonic bonds to morse potentials");
1932 EETYPE("morse", opts->bMorse, yesno_names);
1934 /* Energy group exclusions */
1935 CCTYPE ("ENERGY GROUP EXCLUSIONS");
1936 CTYPE ("Pairs of energy groups for which all non-bonded interactions are excluded");
1937 STYPE ("energygrp-excl", egpexcl, NULL);
1941 CTYPE ("Number of walls, type, atom types, densities and box-z scale factor for Ewald");
1942 ITYPE ("nwall", ir->nwall, 0);
1943 EETYPE("wall-type", ir->wall_type, ewt_names);
1944 RTYPE ("wall-r-linpot", ir->wall_r_linpot, -1);
1945 STYPE ("wall-atomtype", wall_atomtype, NULL);
1946 STYPE ("wall-density", wall_density, NULL);
1947 RTYPE ("wall-ewald-zfac", ir->wall_ewald_zfac, 3);
1950 CCTYPE("COM PULLING");
1951 CTYPE("Pull type: no, umbrella, constraint or constant-force");
1952 EETYPE("pull", ir->ePull, epull_names);
1953 if (ir->ePull != epullNO)
1956 pull_grp = read_pullparams(&ninp, &inp, ir->pull, &opts->pull_start, wi);
1959 /* Enforced rotation */
1960 CCTYPE("ENFORCED ROTATION");
1961 CTYPE("Enforced rotation: No or Yes");
1962 EETYPE("rotation", ir->bRot, yesno_names);
1966 rot_grp = read_rotparams(&ninp, &inp, ir->rot, wi);
1970 CCTYPE("NMR refinement stuff");
1971 CTYPE ("Distance restraints type: No, Simple or Ensemble");
1972 EETYPE("disre", ir->eDisre, edisre_names);
1973 CTYPE ("Force weighting of pairs in one distance restraint: Conservative or Equal");
1974 EETYPE("disre-weighting", ir->eDisreWeighting, edisreweighting_names);
1975 CTYPE ("Use sqrt of the time averaged times the instantaneous violation");
1976 EETYPE("disre-mixed", ir->bDisreMixed, yesno_names);
1977 RTYPE ("disre-fc", ir->dr_fc, 1000.0);
1978 RTYPE ("disre-tau", ir->dr_tau, 0.0);
1979 CTYPE ("Output frequency for pair distances to energy file");
1980 ITYPE ("nstdisreout", ir->nstdisreout, 100);
1981 CTYPE ("Orientation restraints: No or Yes");
1982 EETYPE("orire", opts->bOrire, yesno_names);
1983 CTYPE ("Orientation restraints force constant and tau for time averaging");
1984 RTYPE ("orire-fc", ir->orires_fc, 0.0);
1985 RTYPE ("orire-tau", ir->orires_tau, 0.0);
1986 STYPE ("orire-fitgrp", orirefitgrp, NULL);
1987 CTYPE ("Output frequency for trace(SD) and S to energy file");
1988 ITYPE ("nstorireout", ir->nstorireout, 100);
1990 /* free energy variables */
1991 CCTYPE ("Free energy variables");
1992 EETYPE("free-energy", ir->efep, efep_names);
1993 STYPE ("couple-moltype", couple_moltype, NULL);
1994 EETYPE("couple-lambda0", opts->couple_lam0, couple_lam);
1995 EETYPE("couple-lambda1", opts->couple_lam1, couple_lam);
1996 EETYPE("couple-intramol", opts->bCoupleIntra, yesno_names);
1998 RTYPE ("init-lambda", fep->init_lambda, -1); /* start with -1 so
2000 it was not entered */
2001 ITYPE ("init-lambda-state", fep->init_fep_state, -1);
2002 RTYPE ("delta-lambda", fep->delta_lambda, 0.0);
2003 ITYPE ("nstdhdl", fep->nstdhdl, 50);
2004 STYPE ("fep-lambdas", fep_lambda[efptFEP], NULL);
2005 STYPE ("mass-lambdas", fep_lambda[efptMASS], NULL);
2006 STYPE ("coul-lambdas", fep_lambda[efptCOUL], NULL);
2007 STYPE ("vdw-lambdas", fep_lambda[efptVDW], NULL);
2008 STYPE ("bonded-lambdas", fep_lambda[efptBONDED], NULL);
2009 STYPE ("restraint-lambdas", fep_lambda[efptRESTRAINT], NULL);
2010 STYPE ("temperature-lambdas", fep_lambda[efptTEMPERATURE], NULL);
2011 ITYPE ("calc-lambda-neighbors", fep->lambda_neighbors, 1);
2012 STYPE ("init-lambda-weights", lambda_weights, NULL);
2013 EETYPE("dhdl-print-energy", fep->bPrintEnergy, yesno_names);
2014 RTYPE ("sc-alpha", fep->sc_alpha, 0.0);
2015 ITYPE ("sc-power", fep->sc_power, 1);
2016 RTYPE ("sc-r-power", fep->sc_r_power, 6.0);
2017 RTYPE ("sc-sigma", fep->sc_sigma, 0.3);
2018 EETYPE("sc-coul", fep->bScCoul, yesno_names);
2019 ITYPE ("dh_hist_size", fep->dh_hist_size, 0);
2020 RTYPE ("dh_hist_spacing", fep->dh_hist_spacing, 0.1);
2021 EETYPE("separate-dhdl-file", fep->separate_dhdl_file,
2022 separate_dhdl_file_names);
2023 EETYPE("dhdl-derivatives", fep->dhdl_derivatives, dhdl_derivatives_names);
2024 ITYPE ("dh_hist_size", fep->dh_hist_size, 0);
2025 RTYPE ("dh_hist_spacing", fep->dh_hist_spacing, 0.1);
2027 /* Non-equilibrium MD stuff */
2028 CCTYPE("Non-equilibrium MD stuff");
2029 STYPE ("acc-grps", accgrps, NULL);
2030 STYPE ("accelerate", acc, NULL);
2031 STYPE ("freezegrps", freeze, NULL);
2032 STYPE ("freezedim", frdim, NULL);
2033 RTYPE ("cos-acceleration", ir->cos_accel, 0);
2034 STYPE ("deform", deform, NULL);
2036 /* simulated tempering variables */
2037 CCTYPE("simulated tempering variables");
2038 EETYPE("simulated-tempering", ir->bSimTemp, yesno_names);
2039 EETYPE("simulated-tempering-scaling", ir->simtempvals->eSimTempScale, esimtemp_names);
2040 RTYPE("sim-temp-low", ir->simtempvals->simtemp_low, 300.0);
2041 RTYPE("sim-temp-high", ir->simtempvals->simtemp_high, 300.0);
2043 /* expanded ensemble variables */
2044 if (ir->efep == efepEXPANDED || ir->bSimTemp)
2046 read_expandedparams(&ninp, &inp, expand, wi);
2049 /* Electric fields */
2050 CCTYPE("Electric fields");
2051 CTYPE ("Format is number of terms (int) and for all terms an amplitude (real)");
2052 CTYPE ("and a phase angle (real)");
2053 STYPE ("E-x", efield_x, NULL);
2054 STYPE ("E-xt", efield_xt, NULL);
2055 STYPE ("E-y", efield_y, NULL);
2056 STYPE ("E-yt", efield_yt, NULL);
2057 STYPE ("E-z", efield_z, NULL);
2058 STYPE ("E-zt", efield_zt, NULL);
2060 /* AdResS defined thingies */
2061 CCTYPE ("AdResS parameters");
2062 EETYPE("adress", ir->bAdress, yesno_names);
2065 snew(ir->adress, 1);
2066 read_adressparams(&ninp, &inp, ir->adress, wi);
2069 /* User defined thingies */
2070 CCTYPE ("User defined thingies");
2071 STYPE ("user1-grps", user1, NULL);
2072 STYPE ("user2-grps", user2, NULL);
2073 ITYPE ("userint1", ir->userint1, 0);
2074 ITYPE ("userint2", ir->userint2, 0);
2075 ITYPE ("userint3", ir->userint3, 0);
2076 ITYPE ("userint4", ir->userint4, 0);
2077 RTYPE ("userreal1", ir->userreal1, 0);
2078 RTYPE ("userreal2", ir->userreal2, 0);
2079 RTYPE ("userreal3", ir->userreal3, 0);
2080 RTYPE ("userreal4", ir->userreal4, 0);
2083 write_inpfile(mdparout, ninp, inp, FALSE, wi);
2084 for (i = 0; (i < ninp); i++)
2087 sfree(inp[i].value);
2091 /* Process options if necessary */
2092 for (m = 0; m < 2; m++)
2094 for (i = 0; i < 2*DIM; i++)
2103 if (sscanf(dumstr[m], "%lf", &(dumdub[m][XX])) != 1)
2105 warning_error(wi, "Pressure coupling not enough values (I need 1)");
2107 dumdub[m][YY] = dumdub[m][ZZ] = dumdub[m][XX];
2109 case epctSEMIISOTROPIC:
2110 case epctSURFACETENSION:
2111 if (sscanf(dumstr[m], "%lf%lf",
2112 &(dumdub[m][XX]), &(dumdub[m][ZZ])) != 2)
2114 warning_error(wi, "Pressure coupling not enough values (I need 2)");
2116 dumdub[m][YY] = dumdub[m][XX];
2118 case epctANISOTROPIC:
2119 if (sscanf(dumstr[m], "%lf%lf%lf%lf%lf%lf",
2120 &(dumdub[m][XX]), &(dumdub[m][YY]), &(dumdub[m][ZZ]),
2121 &(dumdub[m][3]), &(dumdub[m][4]), &(dumdub[m][5])) != 6)
2123 warning_error(wi, "Pressure coupling not enough values (I need 6)");
2127 gmx_fatal(FARGS, "Pressure coupling type %s not implemented yet",
2128 epcoupltype_names[ir->epct]);
2132 clear_mat(ir->ref_p);
2133 clear_mat(ir->compress);
2134 for (i = 0; i < DIM; i++)
2136 ir->ref_p[i][i] = dumdub[1][i];
2137 ir->compress[i][i] = dumdub[0][i];
2139 if (ir->epct == epctANISOTROPIC)
2141 ir->ref_p[XX][YY] = dumdub[1][3];
2142 ir->ref_p[XX][ZZ] = dumdub[1][4];
2143 ir->ref_p[YY][ZZ] = dumdub[1][5];
2144 if (ir->ref_p[XX][YY] != 0 && ir->ref_p[XX][ZZ] != 0 && ir->ref_p[YY][ZZ] != 0)
2146 warning(wi, "All off-diagonal reference pressures are non-zero. Are you sure you want to apply a threefold shear stress?\n");
2148 ir->compress[XX][YY] = dumdub[0][3];
2149 ir->compress[XX][ZZ] = dumdub[0][4];
2150 ir->compress[YY][ZZ] = dumdub[0][5];
2151 for (i = 0; i < DIM; i++)
2153 for (m = 0; m < i; m++)
2155 ir->ref_p[i][m] = ir->ref_p[m][i];
2156 ir->compress[i][m] = ir->compress[m][i];
2161 if (ir->comm_mode == ecmNO)
2166 opts->couple_moltype = NULL;
2167 if (strlen(couple_moltype) > 0)
2169 if (ir->efep != efepNO)
2171 opts->couple_moltype = strdup(couple_moltype);
2172 if (opts->couple_lam0 == opts->couple_lam1)
2174 warning(wi, "The lambda=0 and lambda=1 states for coupling are identical");
2176 if (ir->eI == eiMD && (opts->couple_lam0 == ecouplamNONE ||
2177 opts->couple_lam1 == ecouplamNONE))
2179 warning(wi, "For proper sampling of the (nearly) decoupled state, stochastic dynamics should be used");
2184 warning(wi, "Can not couple a molecule with free_energy = no");
2187 /* FREE ENERGY AND EXPANDED ENSEMBLE OPTIONS */
2188 if (ir->efep != efepNO)
2190 if (fep->delta_lambda > 0)
2192 ir->efep = efepSLOWGROWTH;
2198 fep->bPrintEnergy = TRUE;
2199 /* always print out the energy to dhdl if we are doing expanded ensemble, since we need the total energy
2200 if the temperature is changing. */
2203 if ((ir->efep != efepNO) || ir->bSimTemp)
2205 ir->bExpanded = FALSE;
2206 if ((ir->efep == efepEXPANDED) || ir->bSimTemp)
2208 ir->bExpanded = TRUE;
2210 do_fep_params(ir, fep_lambda, lambda_weights);
2211 if (ir->bSimTemp) /* done after fep params */
2213 do_simtemp_params(ir);
2218 ir->fepvals->n_lambda = 0;
2221 /* WALL PARAMETERS */
2223 do_wall_params(ir, wall_atomtype, wall_density, opts);
2225 /* ORIENTATION RESTRAINT PARAMETERS */
2227 if (opts->bOrire && str_nelem(orirefitgrp, MAXPTR, NULL) != 1)
2229 warning_error(wi, "ERROR: Need one orientation restraint fit group\n");
2232 /* DEFORMATION PARAMETERS */
2234 clear_mat(ir->deform);
2235 for (i = 0; i < 6; i++)
2239 m = sscanf(deform, "%lf %lf %lf %lf %lf %lf",
2240 &(dumdub[0][0]), &(dumdub[0][1]), &(dumdub[0][2]),
2241 &(dumdub[0][3]), &(dumdub[0][4]), &(dumdub[0][5]));
2242 for (i = 0; i < 3; i++)
2244 ir->deform[i][i] = dumdub[0][i];
2246 ir->deform[YY][XX] = dumdub[0][3];
2247 ir->deform[ZZ][XX] = dumdub[0][4];
2248 ir->deform[ZZ][YY] = dumdub[0][5];
2249 if (ir->epc != epcNO)
2251 for (i = 0; i < 3; i++)
2253 for (j = 0; j <= i; j++)
2255 if (ir->deform[i][j] != 0 && ir->compress[i][j] != 0)
2257 warning_error(wi, "A box element has deform set and compressibility > 0");
2261 for (i = 0; i < 3; i++)
2263 for (j = 0; j < i; j++)
2265 if (ir->deform[i][j] != 0)
2267 for (m = j; m < DIM; m++)
2269 if (ir->compress[m][j] != 0)
2271 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.");
2272 warning(wi, warn_buf);
2284 static int search_QMstring(char *s, int ng, const char *gn[])
2286 /* same as normal search_string, but this one searches QM strings */
2289 for (i = 0; (i < ng); i++)
2291 if (gmx_strcasecmp(s, gn[i]) == 0)
2297 gmx_fatal(FARGS, "this QM method or basisset (%s) is not implemented\n!", s);
2301 } /* search_QMstring */
2304 int search_string(char *s, int ng, char *gn[])
2308 for (i = 0; (i < ng); i++)
2310 if (gmx_strcasecmp(s, gn[i]) == 0)
2317 "Group %s referenced in the .mdp file was not found in the index file.\n"
2318 "Group names must match either [moleculetype] names or custom index group\n"
2319 "names, in which case you must supply an index file to the '-n' option\n"
2326 static gmx_bool do_numbering(int natoms, gmx_groups_t *groups, int ng, char *ptrs[],
2327 t_blocka *block, char *gnames[],
2328 int gtype, int restnm,
2329 int grptp, gmx_bool bVerbose,
2332 unsigned short *cbuf;
2333 t_grps *grps = &(groups->grps[gtype]);
2334 int i, j, gid, aj, ognr, ntot = 0;
2337 char warn_buf[STRLEN];
2341 fprintf(debug, "Starting numbering %d groups of type %d\n", ng, gtype);
2344 title = gtypes[gtype];
2347 /* Mark all id's as not set */
2348 for (i = 0; (i < natoms); i++)
2353 snew(grps->nm_ind, ng+1); /* +1 for possible rest group */
2354 for (i = 0; (i < ng); i++)
2356 /* Lookup the group name in the block structure */
2357 gid = search_string(ptrs[i], block->nr, gnames);
2358 if ((grptp != egrptpONE) || (i == 0))
2360 grps->nm_ind[grps->nr++] = gid;
2364 fprintf(debug, "Found gid %d for group %s\n", gid, ptrs[i]);
2367 /* Now go over the atoms in the group */
2368 for (j = block->index[gid]; (j < block->index[gid+1]); j++)
2373 /* Range checking */
2374 if ((aj < 0) || (aj >= natoms))
2376 gmx_fatal(FARGS, "Invalid atom number %d in indexfile", aj);
2378 /* Lookup up the old group number */
2382 gmx_fatal(FARGS, "Atom %d in multiple %s groups (%d and %d)",
2383 aj+1, title, ognr+1, i+1);
2387 /* Store the group number in buffer */
2388 if (grptp == egrptpONE)
2401 /* Now check whether we have done all atoms */
2405 if (grptp == egrptpALL)
2407 gmx_fatal(FARGS, "%d atoms are not part of any of the %s groups",
2408 natoms-ntot, title);
2410 else if (grptp == egrptpPART)
2412 sprintf(warn_buf, "%d atoms are not part of any of the %s groups",
2413 natoms-ntot, title);
2414 warning_note(wi, warn_buf);
2416 /* Assign all atoms currently unassigned to a rest group */
2417 for (j = 0; (j < natoms); j++)
2419 if (cbuf[j] == NOGID)
2425 if (grptp != egrptpPART)
2430 "Making dummy/rest group for %s containing %d elements\n",
2431 title, natoms-ntot);
2433 /* Add group name "rest" */
2434 grps->nm_ind[grps->nr] = restnm;
2436 /* Assign the rest name to all atoms not currently assigned to a group */
2437 for (j = 0; (j < natoms); j++)
2439 if (cbuf[j] == NOGID)
2448 if (grps->nr == 1 && (ntot == 0 || ntot == natoms))
2450 /* All atoms are part of one (or no) group, no index required */
2451 groups->ngrpnr[gtype] = 0;
2452 groups->grpnr[gtype] = NULL;
2456 groups->ngrpnr[gtype] = natoms;
2457 snew(groups->grpnr[gtype], natoms);
2458 for (j = 0; (j < natoms); j++)
2460 groups->grpnr[gtype][j] = cbuf[j];
2466 return (bRest && grptp == egrptpPART);
2469 static void calc_nrdf(gmx_mtop_t *mtop, t_inputrec *ir, char **gnames)
2472 gmx_groups_t *groups;
2474 int natoms, ai, aj, i, j, d, g, imin, jmin, nc;
2476 int *nrdf2, *na_vcm, na_tot;
2477 double *nrdf_tc, *nrdf_vcm, nrdf_uc, n_sub = 0;
2478 gmx_mtop_atomloop_all_t aloop;
2480 int mb, mol, ftype, as;
2481 gmx_molblock_t *molb;
2482 gmx_moltype_t *molt;
2485 * First calc 3xnr-atoms for each group
2486 * then subtract half a degree of freedom for each constraint
2488 * Only atoms and nuclei contribute to the degrees of freedom...
2493 groups = &mtop->groups;
2494 natoms = mtop->natoms;
2496 /* Allocate one more for a possible rest group */
2497 /* We need to sum degrees of freedom into doubles,
2498 * since floats give too low nrdf's above 3 million atoms.
2500 snew(nrdf_tc, groups->grps[egcTC].nr+1);
2501 snew(nrdf_vcm, groups->grps[egcVCM].nr+1);
2502 snew(na_vcm, groups->grps[egcVCM].nr+1);
2504 for (i = 0; i < groups->grps[egcTC].nr; i++)
2508 for (i = 0; i < groups->grps[egcVCM].nr+1; i++)
2513 snew(nrdf2, natoms);
2514 aloop = gmx_mtop_atomloop_all_init(mtop);
2515 while (gmx_mtop_atomloop_all_next(aloop, &i, &atom))
2518 if (atom->ptype == eptAtom || atom->ptype == eptNucleus)
2520 g = ggrpnr(groups, egcFREEZE, i);
2521 /* Double count nrdf for particle i */
2522 for (d = 0; d < DIM; d++)
2524 if (opts->nFreeze[g][d] == 0)
2529 nrdf_tc [ggrpnr(groups, egcTC, i)] += 0.5*nrdf2[i];
2530 nrdf_vcm[ggrpnr(groups, egcVCM, i)] += 0.5*nrdf2[i];
2535 for (mb = 0; mb < mtop->nmolblock; mb++)
2537 molb = &mtop->molblock[mb];
2538 molt = &mtop->moltype[molb->type];
2539 atom = molt->atoms.atom;
2540 for (mol = 0; mol < molb->nmol; mol++)
2542 for (ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++)
2544 ia = molt->ilist[ftype].iatoms;
2545 for (i = 0; i < molt->ilist[ftype].nr; )
2547 /* Subtract degrees of freedom for the constraints,
2548 * if the particles still have degrees of freedom left.
2549 * If one of the particles is a vsite or a shell, then all
2550 * constraint motion will go there, but since they do not
2551 * contribute to the constraints the degrees of freedom do not
2556 if (((atom[ia[1]].ptype == eptNucleus) ||
2557 (atom[ia[1]].ptype == eptAtom)) &&
2558 ((atom[ia[2]].ptype == eptNucleus) ||
2559 (atom[ia[2]].ptype == eptAtom)))
2577 imin = min(imin, nrdf2[ai]);
2578 jmin = min(jmin, nrdf2[aj]);
2581 nrdf_tc [ggrpnr(groups, egcTC, ai)] -= 0.5*imin;
2582 nrdf_tc [ggrpnr(groups, egcTC, aj)] -= 0.5*jmin;
2583 nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5*imin;
2584 nrdf_vcm[ggrpnr(groups, egcVCM, aj)] -= 0.5*jmin;
2586 ia += interaction_function[ftype].nratoms+1;
2587 i += interaction_function[ftype].nratoms+1;
2590 ia = molt->ilist[F_SETTLE].iatoms;
2591 for (i = 0; i < molt->ilist[F_SETTLE].nr; )
2593 /* Subtract 1 dof from every atom in the SETTLE */
2594 for (j = 0; j < 3; j++)
2597 imin = min(2, nrdf2[ai]);
2599 nrdf_tc [ggrpnr(groups, egcTC, ai)] -= 0.5*imin;
2600 nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5*imin;
2605 as += molt->atoms.nr;
2609 if (ir->ePull == epullCONSTRAINT)
2611 /* Correct nrdf for the COM constraints.
2612 * We correct using the TC and VCM group of the first atom
2613 * in the reference and pull group. If atoms in one pull group
2614 * belong to different TC or VCM groups it is anyhow difficult
2615 * to determine the optimal nrdf assignment.
2618 if (pull->eGeom == epullgPOS)
2621 for (i = 0; i < DIM; i++)
2633 for (i = 0; i < pull->ngrp; i++)
2636 if (pull->grp[0].nat > 0)
2638 /* Subtract 1/2 dof from the reference group */
2639 ai = pull->grp[0].ind[0];
2640 if (nrdf_tc[ggrpnr(groups, egcTC, ai)] > 1)
2642 nrdf_tc [ggrpnr(groups, egcTC, ai)] -= 0.5;
2643 nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5;
2647 /* Subtract 1/2 dof from the pulled group */
2648 ai = pull->grp[1+i].ind[0];
2649 nrdf_tc [ggrpnr(groups, egcTC, ai)] -= 0.5*imin;
2650 nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5*imin;
2651 if (nrdf_tc[ggrpnr(groups, egcTC, ai)] < 0)
2653 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)]]);
2658 if (ir->nstcomm != 0)
2660 /* Subtract 3 from the number of degrees of freedom in each vcm group
2661 * when com translation is removed and 6 when rotation is removed
2664 switch (ir->comm_mode)
2667 n_sub = ndof_com(ir);
2674 gmx_incons("Checking comm_mode");
2677 for (i = 0; i < groups->grps[egcTC].nr; i++)
2679 /* Count the number of atoms of TC group i for every VCM group */
2680 for (j = 0; j < groups->grps[egcVCM].nr+1; j++)
2685 for (ai = 0; ai < natoms; ai++)
2687 if (ggrpnr(groups, egcTC, ai) == i)
2689 na_vcm[ggrpnr(groups, egcVCM, ai)]++;
2693 /* Correct for VCM removal according to the fraction of each VCM
2694 * group present in this TC group.
2696 nrdf_uc = nrdf_tc[i];
2699 fprintf(debug, "T-group[%d] nrdf_uc = %g, n_sub = %g\n",
2703 for (j = 0; j < groups->grps[egcVCM].nr+1; j++)
2705 if (nrdf_vcm[j] > n_sub)
2707 nrdf_tc[i] += nrdf_uc*((double)na_vcm[j]/(double)na_tot)*
2708 (nrdf_vcm[j] - n_sub)/nrdf_vcm[j];
2712 fprintf(debug, " nrdf_vcm[%d] = %g, nrdf = %g\n",
2713 j, nrdf_vcm[j], nrdf_tc[i]);
2718 for (i = 0; (i < groups->grps[egcTC].nr); i++)
2720 opts->nrdf[i] = nrdf_tc[i];
2721 if (opts->nrdf[i] < 0)
2726 "Number of degrees of freedom in T-Coupling group %s is %.2f\n",
2727 gnames[groups->grps[egcTC].nm_ind[i]], opts->nrdf[i]);
2736 static void decode_cos(char *s, t_cosines *cosine, gmx_bool bTime)
2739 char format[STRLEN], f1[STRLEN];
2751 sscanf(t, "%d", &(cosine->n));
2758 snew(cosine->a, cosine->n);
2759 snew(cosine->phi, cosine->n);
2761 sprintf(format, "%%*d");
2762 for (i = 0; (i < cosine->n); i++)
2765 strcat(f1, "%lf%lf");
2766 if (sscanf(t, f1, &a, &phi) < 2)
2768 gmx_fatal(FARGS, "Invalid input for electric field shift: '%s'", t);
2771 cosine->phi[i] = phi;
2772 strcat(format, "%*lf%*lf");
2779 static gmx_bool do_egp_flag(t_inputrec *ir, gmx_groups_t *groups,
2780 const char *option, const char *val, int flag)
2782 /* The maximum number of energy group pairs would be MAXPTR*(MAXPTR+1)/2.
2783 * But since this is much larger than STRLEN, such a line can not be parsed.
2784 * The real maximum is the number of names that fit in a string: STRLEN/2.
2786 #define EGP_MAX (STRLEN/2)
2787 int nelem, i, j, k, nr;
2788 char *names[EGP_MAX];
2792 gnames = groups->grpname;
2794 nelem = str_nelem(val, EGP_MAX, names);
2797 gmx_fatal(FARGS, "The number of groups for %s is odd", option);
2799 nr = groups->grps[egcENER].nr;
2801 for (i = 0; i < nelem/2; i++)
2805 gmx_strcasecmp(names[2*i], *(gnames[groups->grps[egcENER].nm_ind[j]])))
2811 gmx_fatal(FARGS, "%s in %s is not an energy group\n",
2812 names[2*i], option);
2816 gmx_strcasecmp(names[2*i+1], *(gnames[groups->grps[egcENER].nm_ind[k]])))
2822 gmx_fatal(FARGS, "%s in %s is not an energy group\n",
2823 names[2*i+1], option);
2825 if ((j < nr) && (k < nr))
2827 ir->opts.egp_flags[nr*j+k] |= flag;
2828 ir->opts.egp_flags[nr*k+j] |= flag;
2836 void do_index(const char* mdparin, const char *ndx,
2839 t_inputrec *ir, rvec *v,
2843 gmx_groups_t *groups;
2847 char warnbuf[STRLEN], **gnames;
2848 int nr, ntcg, ntau_t, nref_t, nacc, nofg, nSA, nSA_points, nSA_time, nSA_temp;
2851 int nacg, nfreeze, nfrdim, nenergy, nvcm, nuser;
2852 char *ptr1[MAXPTR], *ptr2[MAXPTR], *ptr3[MAXPTR];
2853 int i, j, k, restnm;
2855 gmx_bool bExcl, bTable, bSetTCpar, bAnneal, bRest;
2856 int nQMmethod, nQMbasis, nQMcharge, nQMmult, nbSH, nCASorb, nCASelec,
2857 nSAon, nSAoff, nSAsteps, nQMg, nbOPT, nbTS;
2858 char warn_buf[STRLEN];
2862 fprintf(stderr, "processing index file...\n");
2868 snew(grps->index, 1);
2870 atoms_all = gmx_mtop_global_atoms(mtop);
2871 analyse(&atoms_all, grps, &gnames, FALSE, TRUE);
2872 free_t_atoms(&atoms_all, FALSE);
2876 grps = init_index(ndx, &gnames);
2879 groups = &mtop->groups;
2880 natoms = mtop->natoms;
2881 symtab = &mtop->symtab;
2883 snew(groups->grpname, grps->nr+1);
2885 for (i = 0; (i < grps->nr); i++)
2887 groups->grpname[i] = put_symtab(symtab, gnames[i]);
2889 groups->grpname[i] = put_symtab(symtab, "rest");
2891 srenew(gnames, grps->nr+1);
2892 gnames[restnm] = *(groups->grpname[i]);
2893 groups->ngrpname = grps->nr+1;
2895 set_warning_line(wi, mdparin, -1);
2897 ntau_t = str_nelem(tau_t, MAXPTR, ptr1);
2898 nref_t = str_nelem(ref_t, MAXPTR, ptr2);
2899 ntcg = str_nelem(tcgrps, MAXPTR, ptr3);
2900 if ((ntau_t != ntcg) || (nref_t != ntcg))
2902 gmx_fatal(FARGS, "Invalid T coupling input: %d groups, %d ref-t values and "
2903 "%d tau-t values", ntcg, nref_t, ntau_t);
2906 bSetTCpar = (ir->etc || EI_SD(ir->eI) || ir->eI == eiBD || EI_TPI(ir->eI));
2907 do_numbering(natoms, groups, ntcg, ptr3, grps, gnames, egcTC,
2908 restnm, bSetTCpar ? egrptpALL : egrptpALL_GENREST, bVerbose, wi);
2909 nr = groups->grps[egcTC].nr;
2911 snew(ir->opts.nrdf, nr);
2912 snew(ir->opts.tau_t, nr);
2913 snew(ir->opts.ref_t, nr);
2914 if (ir->eI == eiBD && ir->bd_fric == 0)
2916 fprintf(stderr, "bd-fric=0, so tau-t will be used as the inverse friction constant(s)\n");
2923 gmx_fatal(FARGS, "Not enough ref-t and tau-t values!");
2927 for (i = 0; (i < nr); i++)
2929 ir->opts.tau_t[i] = strtod(ptr1[i], NULL);
2930 if ((ir->eI == eiBD || ir->eI == eiSD2) && ir->opts.tau_t[i] <= 0)
2932 sprintf(warn_buf, "With integrator %s tau-t should be larger than 0", ei_names[ir->eI]);
2933 warning_error(wi, warn_buf);
2936 if (ir->etc != etcVRESCALE && ir->opts.tau_t[i] == 0)
2938 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.");
2941 if (ir->opts.tau_t[i] >= 0)
2943 tau_min = min(tau_min, ir->opts.tau_t[i]);
2946 if (ir->etc != etcNO && ir->nsttcouple == -1)
2948 ir->nsttcouple = ir_optimal_nsttcouple(ir);
2953 if ((ir->etc == etcNOSEHOOVER) && (ir->epc == epcBERENDSEN))
2955 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");
2957 if ((ir->epc == epcMTTK) && (ir->etc > etcNO))
2959 if (ir->nstpcouple != ir->nsttcouple)
2961 int mincouple = min(ir->nstpcouple, ir->nsttcouple);
2962 ir->nstpcouple = ir->nsttcouple = mincouple;
2963 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);
2964 warning_note(wi, warn_buf);
2968 /* velocity verlet with averaged kinetic energy KE = 0.5*(v(t+1/2) - v(t-1/2)) is implemented
2969 primarily for testing purposes, and does not work with temperature coupling other than 1 */
2971 if (ETC_ANDERSEN(ir->etc))
2973 if (ir->nsttcouple != 1)
2976 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");
2977 warning_note(wi, warn_buf);
2980 nstcmin = tcouple_min_integration_steps(ir->etc);
2983 if (tau_min/(ir->delta_t*ir->nsttcouple) < nstcmin)
2985 sprintf(warn_buf, "For proper integration of the %s thermostat, tau-t (%g) should be at least %d times larger than nsttcouple*dt (%g)",
2986 ETCOUPLTYPE(ir->etc),
2988 ir->nsttcouple*ir->delta_t);
2989 warning(wi, warn_buf);
2992 for (i = 0; (i < nr); i++)
2994 ir->opts.ref_t[i] = strtod(ptr2[i], NULL);
2995 if (ir->opts.ref_t[i] < 0)
2997 gmx_fatal(FARGS, "ref-t for group %d negative", i);
3000 /* set the lambda mc temperature to the md integrator temperature (which should be defined
3001 if we are in this conditional) if mc_temp is negative */
3002 if (ir->expandedvals->mc_temp < 0)
3004 ir->expandedvals->mc_temp = ir->opts.ref_t[0]; /*for now, set to the first reft */
3008 /* Simulated annealing for each group. There are nr groups */
3009 nSA = str_nelem(anneal, MAXPTR, ptr1);
3010 if (nSA == 1 && (ptr1[0][0] == 'n' || ptr1[0][0] == 'N'))
3014 if (nSA > 0 && nSA != nr)
3016 gmx_fatal(FARGS, "Not enough annealing values: %d (for %d groups)\n", nSA, nr);
3020 snew(ir->opts.annealing, nr);
3021 snew(ir->opts.anneal_npoints, nr);
3022 snew(ir->opts.anneal_time, nr);
3023 snew(ir->opts.anneal_temp, nr);
3024 for (i = 0; i < nr; i++)
3026 ir->opts.annealing[i] = eannNO;
3027 ir->opts.anneal_npoints[i] = 0;
3028 ir->opts.anneal_time[i] = NULL;
3029 ir->opts.anneal_temp[i] = NULL;
3034 for (i = 0; i < nr; i++)
3036 if (ptr1[i][0] == 'n' || ptr1[i][0] == 'N')
3038 ir->opts.annealing[i] = eannNO;
3040 else if (ptr1[i][0] == 's' || ptr1[i][0] == 'S')
3042 ir->opts.annealing[i] = eannSINGLE;
3045 else if (ptr1[i][0] == 'p' || ptr1[i][0] == 'P')
3047 ir->opts.annealing[i] = eannPERIODIC;
3053 /* Read the other fields too */
3054 nSA_points = str_nelem(anneal_npoints, MAXPTR, ptr1);
3055 if (nSA_points != nSA)
3057 gmx_fatal(FARGS, "Found %d annealing-npoints values for %d groups\n", nSA_points, nSA);
3059 for (k = 0, i = 0; i < nr; i++)
3061 ir->opts.anneal_npoints[i] = strtol(ptr1[i], NULL, 10);
3062 if (ir->opts.anneal_npoints[i] == 1)
3064 gmx_fatal(FARGS, "Please specify at least a start and an end point for annealing\n");
3066 snew(ir->opts.anneal_time[i], ir->opts.anneal_npoints[i]);
3067 snew(ir->opts.anneal_temp[i], ir->opts.anneal_npoints[i]);
3068 k += ir->opts.anneal_npoints[i];
3071 nSA_time = str_nelem(anneal_time, MAXPTR, ptr1);
3074 gmx_fatal(FARGS, "Found %d annealing-time values, wanter %d\n", nSA_time, k);
3076 nSA_temp = str_nelem(anneal_temp, MAXPTR, ptr2);
3079 gmx_fatal(FARGS, "Found %d annealing-temp values, wanted %d\n", nSA_temp, k);
3082 for (i = 0, k = 0; i < nr; i++)
3085 for (j = 0; j < ir->opts.anneal_npoints[i]; j++)
3087 ir->opts.anneal_time[i][j] = strtod(ptr1[k], NULL);
3088 ir->opts.anneal_temp[i][j] = strtod(ptr2[k], NULL);
3091 if (ir->opts.anneal_time[i][0] > (ir->init_t+GMX_REAL_EPS))
3093 gmx_fatal(FARGS, "First time point for annealing > init_t.\n");
3099 if (ir->opts.anneal_time[i][j] < ir->opts.anneal_time[i][j-1])
3101 gmx_fatal(FARGS, "Annealing timepoints out of order: t=%f comes after t=%f\n",
3102 ir->opts.anneal_time[i][j], ir->opts.anneal_time[i][j-1]);
3105 if (ir->opts.anneal_temp[i][j] < 0)
3107 gmx_fatal(FARGS, "Found negative temperature in annealing: %f\n", ir->opts.anneal_temp[i][j]);
3112 /* Print out some summary information, to make sure we got it right */
3113 for (i = 0, k = 0; i < nr; i++)
3115 if (ir->opts.annealing[i] != eannNO)
3117 j = groups->grps[egcTC].nm_ind[i];
3118 fprintf(stderr, "Simulated annealing for group %s: %s, %d timepoints\n",
3119 *(groups->grpname[j]), eann_names[ir->opts.annealing[i]],
3120 ir->opts.anneal_npoints[i]);
3121 fprintf(stderr, "Time (ps) Temperature (K)\n");
3122 /* All terms except the last one */
3123 for (j = 0; j < (ir->opts.anneal_npoints[i]-1); j++)
3125 fprintf(stderr, "%9.1f %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3128 /* Finally the last one */
3129 j = ir->opts.anneal_npoints[i]-1;
3130 if (ir->opts.annealing[i] == eannSINGLE)
3132 fprintf(stderr, "%9.1f- %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3136 fprintf(stderr, "%9.1f %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3137 if (fabs(ir->opts.anneal_temp[i][j]-ir->opts.anneal_temp[i][0]) > GMX_REAL_EPS)
3139 warning_note(wi, "There is a temperature jump when your annealing loops back.\n");
3148 if (ir->ePull != epullNO)
3150 make_pull_groups(ir->pull, pull_grp, grps, gnames);
3155 make_rotation_groups(ir->rot, rot_grp, grps, gnames);
3158 nacc = str_nelem(acc, MAXPTR, ptr1);
3159 nacg = str_nelem(accgrps, MAXPTR, ptr2);
3160 if (nacg*DIM != nacc)
3162 gmx_fatal(FARGS, "Invalid Acceleration input: %d groups and %d acc. values",
3165 do_numbering(natoms, groups, nacg, ptr2, grps, gnames, egcACC,
3166 restnm, egrptpALL_GENREST, bVerbose, wi);
3167 nr = groups->grps[egcACC].nr;
3168 snew(ir->opts.acc, nr);
3169 ir->opts.ngacc = nr;
3171 for (i = k = 0; (i < nacg); i++)
3173 for (j = 0; (j < DIM); j++, k++)
3175 ir->opts.acc[i][j] = strtod(ptr1[k], NULL);
3178 for (; (i < nr); i++)
3180 for (j = 0; (j < DIM); j++)
3182 ir->opts.acc[i][j] = 0;
3186 nfrdim = str_nelem(frdim, MAXPTR, ptr1);
3187 nfreeze = str_nelem(freeze, MAXPTR, ptr2);
3188 if (nfrdim != DIM*nfreeze)
3190 gmx_fatal(FARGS, "Invalid Freezing input: %d groups and %d freeze values",
3193 do_numbering(natoms, groups, nfreeze, ptr2, grps, gnames, egcFREEZE,
3194 restnm, egrptpALL_GENREST, bVerbose, wi);
3195 nr = groups->grps[egcFREEZE].nr;
3196 ir->opts.ngfrz = nr;
3197 snew(ir->opts.nFreeze, nr);
3198 for (i = k = 0; (i < nfreeze); i++)
3200 for (j = 0; (j < DIM); j++, k++)
3202 ir->opts.nFreeze[i][j] = (gmx_strncasecmp(ptr1[k], "Y", 1) == 0);
3203 if (!ir->opts.nFreeze[i][j])
3205 if (gmx_strncasecmp(ptr1[k], "N", 1) != 0)
3207 sprintf(warnbuf, "Please use Y(ES) or N(O) for freezedim only "
3208 "(not %s)", ptr1[k]);
3209 warning(wi, warn_buf);
3214 for (; (i < nr); i++)
3216 for (j = 0; (j < DIM); j++)
3218 ir->opts.nFreeze[i][j] = 0;
3222 nenergy = str_nelem(energy, MAXPTR, ptr1);
3223 do_numbering(natoms, groups, nenergy, ptr1, grps, gnames, egcENER,
3224 restnm, egrptpALL_GENREST, bVerbose, wi);
3225 add_wall_energrps(groups, ir->nwall, symtab);
3226 ir->opts.ngener = groups->grps[egcENER].nr;
3227 nvcm = str_nelem(vcm, MAXPTR, ptr1);
3229 do_numbering(natoms, groups, nvcm, ptr1, grps, gnames, egcVCM,
3230 restnm, nvcm == 0 ? egrptpALL_GENREST : egrptpPART, bVerbose, wi);
3233 warning(wi, "Some atoms are not part of any center of mass motion removal group.\n"
3234 "This may lead to artifacts.\n"
3235 "In most cases one should use one group for the whole system.");
3238 /* Now we have filled the freeze struct, so we can calculate NRDF */
3239 calc_nrdf(mtop, ir, gnames);
3245 /* Must check per group! */
3246 for (i = 0; (i < ir->opts.ngtc); i++)
3248 ntot += ir->opts.nrdf[i];
3250 if (ntot != (DIM*natoms))
3252 fac = sqrt(ntot/(DIM*natoms));
3255 fprintf(stderr, "Scaling velocities by a factor of %.3f to account for constraints\n"
3256 "and removal of center of mass motion\n", fac);
3258 for (i = 0; (i < natoms); i++)
3260 svmul(fac, v[i], v[i]);
3265 nuser = str_nelem(user1, MAXPTR, ptr1);
3266 do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcUser1,
3267 restnm, egrptpALL_GENREST, bVerbose, wi);
3268 nuser = str_nelem(user2, MAXPTR, ptr1);
3269 do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcUser2,
3270 restnm, egrptpALL_GENREST, bVerbose, wi);
3271 nuser = str_nelem(xtc_grps, MAXPTR, ptr1);
3272 do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcXTC,
3273 restnm, egrptpONE, bVerbose, wi);
3274 nofg = str_nelem(orirefitgrp, MAXPTR, ptr1);
3275 do_numbering(natoms, groups, nofg, ptr1, grps, gnames, egcORFIT,
3276 restnm, egrptpALL_GENREST, bVerbose, wi);
3278 /* QMMM input processing */
3279 nQMg = str_nelem(QMMM, MAXPTR, ptr1);
3280 nQMmethod = str_nelem(QMmethod, MAXPTR, ptr2);
3281 nQMbasis = str_nelem(QMbasis, MAXPTR, ptr3);
3282 if ((nQMmethod != nQMg) || (nQMbasis != nQMg))
3284 gmx_fatal(FARGS, "Invalid QMMM input: %d groups %d basissets"
3285 " and %d methods\n", nQMg, nQMbasis, nQMmethod);
3287 /* group rest, if any, is always MM! */
3288 do_numbering(natoms, groups, nQMg, ptr1, grps, gnames, egcQMMM,
3289 restnm, egrptpALL_GENREST, bVerbose, wi);
3290 nr = nQMg; /*atoms->grps[egcQMMM].nr;*/
3291 ir->opts.ngQM = nQMg;
3292 snew(ir->opts.QMmethod, nr);
3293 snew(ir->opts.QMbasis, nr);
3294 for (i = 0; i < nr; i++)
3296 /* input consists of strings: RHF CASSCF PM3 .. These need to be
3297 * converted to the corresponding enum in names.c
3299 ir->opts.QMmethod[i] = search_QMstring(ptr2[i], eQMmethodNR,
3301 ir->opts.QMbasis[i] = search_QMstring(ptr3[i], eQMbasisNR,
3305 nQMmult = str_nelem(QMmult, MAXPTR, ptr1);
3306 nQMcharge = str_nelem(QMcharge, MAXPTR, ptr2);
3307 nbSH = str_nelem(bSH, MAXPTR, ptr3);
3308 snew(ir->opts.QMmult, nr);
3309 snew(ir->opts.QMcharge, nr);
3310 snew(ir->opts.bSH, nr);
3312 for (i = 0; i < nr; i++)
3314 ir->opts.QMmult[i] = strtol(ptr1[i], NULL, 10);
3315 ir->opts.QMcharge[i] = strtol(ptr2[i], NULL, 10);
3316 ir->opts.bSH[i] = (gmx_strncasecmp(ptr3[i], "Y", 1) == 0);
3319 nCASelec = str_nelem(CASelectrons, MAXPTR, ptr1);
3320 nCASorb = str_nelem(CASorbitals, MAXPTR, ptr2);
3321 snew(ir->opts.CASelectrons, nr);
3322 snew(ir->opts.CASorbitals, nr);
3323 for (i = 0; i < nr; i++)
3325 ir->opts.CASelectrons[i] = strtol(ptr1[i], NULL, 10);
3326 ir->opts.CASorbitals[i] = strtol(ptr2[i], NULL, 10);
3328 /* special optimization options */
3330 nbOPT = str_nelem(bOPT, MAXPTR, ptr1);
3331 nbTS = str_nelem(bTS, MAXPTR, ptr2);
3332 snew(ir->opts.bOPT, nr);
3333 snew(ir->opts.bTS, nr);
3334 for (i = 0; i < nr; i++)
3336 ir->opts.bOPT[i] = (gmx_strncasecmp(ptr1[i], "Y", 1) == 0);
3337 ir->opts.bTS[i] = (gmx_strncasecmp(ptr2[i], "Y", 1) == 0);
3339 nSAon = str_nelem(SAon, MAXPTR, ptr1);
3340 nSAoff = str_nelem(SAoff, MAXPTR, ptr2);
3341 nSAsteps = str_nelem(SAsteps, MAXPTR, ptr3);
3342 snew(ir->opts.SAon, nr);
3343 snew(ir->opts.SAoff, nr);
3344 snew(ir->opts.SAsteps, nr);
3346 for (i = 0; i < nr; i++)
3348 ir->opts.SAon[i] = strtod(ptr1[i], NULL);
3349 ir->opts.SAoff[i] = strtod(ptr2[i], NULL);
3350 ir->opts.SAsteps[i] = strtol(ptr3[i], NULL, 10);
3352 /* end of QMMM input */
3356 for (i = 0; (i < egcNR); i++)
3358 fprintf(stderr, "%-16s has %d element(s):", gtypes[i], groups->grps[i].nr);
3359 for (j = 0; (j < groups->grps[i].nr); j++)
3361 fprintf(stderr, " %s", *(groups->grpname[groups->grps[i].nm_ind[j]]));
3363 fprintf(stderr, "\n");
3367 nr = groups->grps[egcENER].nr;
3368 snew(ir->opts.egp_flags, nr*nr);
3370 bExcl = do_egp_flag(ir, groups, "energygrp-excl", egpexcl, EGP_EXCL);
3371 if (bExcl && ir->cutoff_scheme == ecutsVERLET)
3373 warning_error(wi, "Energy group exclusions are not (yet) implemented for the Verlet scheme");
3375 if (bExcl && EEL_FULL(ir->coulombtype))
3377 warning(wi, "Can not exclude the lattice Coulomb energy between energy groups");
3380 bTable = do_egp_flag(ir, groups, "energygrp-table", egptable, EGP_TABLE);
3381 if (bTable && !(ir->vdwtype == evdwUSER) &&
3382 !(ir->coulombtype == eelUSER) && !(ir->coulombtype == eelPMEUSER) &&
3383 !(ir->coulombtype == eelPMEUSERSWITCH))
3385 gmx_fatal(FARGS, "Can only have energy group pair tables in combination with user tables for VdW and/or Coulomb");
3388 decode_cos(efield_x, &(ir->ex[XX]), FALSE);
3389 decode_cos(efield_xt, &(ir->et[XX]), TRUE);
3390 decode_cos(efield_y, &(ir->ex[YY]), FALSE);
3391 decode_cos(efield_yt, &(ir->et[YY]), TRUE);
3392 decode_cos(efield_z, &(ir->ex[ZZ]), FALSE);
3393 decode_cos(efield_zt, &(ir->et[ZZ]), TRUE);
3397 do_adress_index(ir->adress, groups, gnames, &(ir->opts), wi);
3400 for (i = 0; (i < grps->nr); i++)
3412 static void check_disre(gmx_mtop_t *mtop)
3414 gmx_ffparams_t *ffparams;
3415 t_functype *functype;
3417 int i, ndouble, ftype;
3418 int label, old_label;
3420 if (gmx_mtop_ftype_count(mtop, F_DISRES) > 0)
3422 ffparams = &mtop->ffparams;
3423 functype = ffparams->functype;
3424 ip = ffparams->iparams;
3427 for (i = 0; i < ffparams->ntypes; i++)
3429 ftype = functype[i];
3430 if (ftype == F_DISRES)
3432 label = ip[i].disres.label;
3433 if (label == old_label)
3435 fprintf(stderr, "Distance restraint index %d occurs twice\n", label);
3443 gmx_fatal(FARGS, "Found %d double distance restraint indices,\n"
3444 "probably the parameters for multiple pairs in one restraint "
3445 "are not identical\n", ndouble);
3450 static gmx_bool absolute_reference(t_inputrec *ir, gmx_mtop_t *sys,
3451 gmx_bool posres_only,
3455 gmx_mtop_ilistloop_t iloop;
3465 for (d = 0; d < DIM; d++)
3467 AbsRef[d] = (d < ndof_com(ir) ? 0 : 1);
3469 /* Check for freeze groups */
3470 for (g = 0; g < ir->opts.ngfrz; g++)
3472 for (d = 0; d < DIM; d++)
3474 if (ir->opts.nFreeze[g][d] != 0)
3482 /* Check for position restraints */
3483 iloop = gmx_mtop_ilistloop_init(sys);
3484 while (gmx_mtop_ilistloop_next(iloop, &ilist, &nmol))
3487 (AbsRef[XX] == 0 || AbsRef[YY] == 0 || AbsRef[ZZ] == 0))
3489 for (i = 0; i < ilist[F_POSRES].nr; i += 2)
3491 pr = &sys->ffparams.iparams[ilist[F_POSRES].iatoms[i]];
3492 for (d = 0; d < DIM; d++)
3494 if (pr->posres.fcA[d] != 0)
3503 return (AbsRef[XX] != 0 && AbsRef[YY] != 0 && AbsRef[ZZ] != 0);
3506 void triple_check(const char *mdparin, t_inputrec *ir, gmx_mtop_t *sys,
3510 int i, m, g, nmol, npct;
3511 gmx_bool bCharge, bAcc;
3512 real gdt_max, *mgrp, mt;
3514 gmx_mtop_atomloop_block_t aloopb;
3515 gmx_mtop_atomloop_all_t aloop;
3518 char warn_buf[STRLEN];
3520 set_warning_line(wi, mdparin, -1);
3522 if (EI_DYNAMICS(ir->eI) && !EI_SD(ir->eI) && ir->eI != eiBD &&
3523 ir->comm_mode == ecmNO &&
3524 !(absolute_reference(ir, sys, FALSE, AbsRef) || ir->nsteps <= 10) &&
3525 !ETC_ANDERSEN(ir->etc))
3527 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");
3530 /* Check for pressure coupling with absolute position restraints */
3531 if (ir->epc != epcNO && ir->refcoord_scaling == erscNO)
3533 absolute_reference(ir, sys, TRUE, AbsRef);
3535 for (m = 0; m < DIM; m++)
3537 if (AbsRef[m] && norm2(ir->compress[m]) > 0)
3539 warning(wi, "You are using pressure coupling with absolute position restraints, this will give artifacts. Use the refcoord_scaling option.");
3547 aloopb = gmx_mtop_atomloop_block_init(sys);
3548 while (gmx_mtop_atomloop_block_next(aloopb, &atom, &nmol))
3550 if (atom->q != 0 || atom->qB != 0)
3558 if (EEL_FULL(ir->coulombtype))
3561 "You are using full electrostatics treatment %s for a system without charges.\n"
3562 "This costs a lot of performance for just processing zeros, consider using %s instead.\n",
3563 EELTYPE(ir->coulombtype), EELTYPE(eelCUT));
3564 warning(wi, err_buf);
3569 if (ir->coulombtype == eelCUT && ir->rcoulomb > 0 && !ir->implicit_solvent)
3572 "You are using a plain Coulomb cut-off, which might produce artifacts.\n"
3573 "You might want to consider using %s electrostatics.\n",
3575 warning_note(wi, err_buf);
3579 /* Generalized reaction field */
3580 if (ir->opts.ngtc == 0)
3582 sprintf(err_buf, "No temperature coupling while using coulombtype %s",
3584 CHECK(ir->coulombtype == eelGRF);
3588 sprintf(err_buf, "When using coulombtype = %s"
3589 " ref-t for temperature coupling should be > 0",
3591 CHECK((ir->coulombtype == eelGRF) && (ir->opts.ref_t[0] <= 0));
3594 if (ir->eI == eiSD1 &&
3595 (gmx_mtop_ftype_count(sys, F_CONSTR) > 0 ||
3596 gmx_mtop_ftype_count(sys, F_SETTLE) > 0))
3598 sprintf(warn_buf, "With constraints integrator %s is less accurate, consider using %s instead", ei_names[ir->eI], ei_names[eiSD2]);
3599 warning_note(wi, warn_buf);
3603 for (i = 0; (i < sys->groups.grps[egcACC].nr); i++)
3605 for (m = 0; (m < DIM); m++)
3607 if (fabs(ir->opts.acc[i][m]) > 1e-6)
3616 snew(mgrp, sys->groups.grps[egcACC].nr);
3617 aloop = gmx_mtop_atomloop_all_init(sys);
3618 while (gmx_mtop_atomloop_all_next(aloop, &i, &atom))
3620 mgrp[ggrpnr(&sys->groups, egcACC, i)] += atom->m;
3623 for (i = 0; (i < sys->groups.grps[egcACC].nr); i++)
3625 for (m = 0; (m < DIM); m++)
3627 acc[m] += ir->opts.acc[i][m]*mgrp[i];
3631 for (m = 0; (m < DIM); m++)
3633 if (fabs(acc[m]) > 1e-6)
3635 const char *dim[DIM] = { "X", "Y", "Z" };
3637 "Net Acceleration in %s direction, will %s be corrected\n",
3638 dim[m], ir->nstcomm != 0 ? "" : "not");
3639 if (ir->nstcomm != 0 && m < ndof_com(ir))
3642 for (i = 0; (i < sys->groups.grps[egcACC].nr); i++)
3644 ir->opts.acc[i][m] -= acc[m];
3652 if (ir->efep != efepNO && ir->fepvals->sc_alpha != 0 &&
3653 !gmx_within_tol(sys->ffparams.reppow, 12.0, 10*GMX_DOUBLE_EPS))
3655 gmx_fatal(FARGS, "Soft-core interactions are only supported with VdW repulsion power 12");
3658 if (ir->ePull != epullNO)
3660 if (ir->pull->grp[0].nat == 0)
3662 absolute_reference(ir, sys, FALSE, AbsRef);
3663 for (m = 0; m < DIM; m++)
3665 if (ir->pull->dim[m] && !AbsRef[m])
3667 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.");
3673 if (ir->pull->eGeom == epullgDIRPBC)
3675 for (i = 0; i < 3; i++)
3677 for (m = 0; m <= i; m++)
3679 if ((ir->epc != epcNO && ir->compress[i][m] != 0) ||
3680 ir->deform[i][m] != 0)
3682 for (g = 1; g < ir->pull->ngrp; g++)
3684 if (ir->pull->grp[g].vec[m] != 0)
3686 gmx_fatal(FARGS, "Can not have dynamic box while using pull geometry '%s' (dim %c)", EPULLGEOM(ir->pull->eGeom), 'x'+m);
3698 void double_check(t_inputrec *ir, matrix box, gmx_bool bConstr, warninp_t wi)
3702 char warn_buf[STRLEN];
3705 ptr = check_box(ir->ePBC, box);
3708 warning_error(wi, ptr);
3711 if (bConstr && ir->eConstrAlg == econtSHAKE)
3713 if (ir->shake_tol <= 0.0)
3715 sprintf(warn_buf, "ERROR: shake-tol must be > 0 instead of %g\n",
3717 warning_error(wi, warn_buf);
3720 if (IR_TWINRANGE(*ir) && ir->nstlist > 1)
3722 sprintf(warn_buf, "With twin-range cut-off's and SHAKE the virial and the pressure are incorrect.");
3723 if (ir->epc == epcNO)
3725 warning(wi, warn_buf);
3729 warning_error(wi, warn_buf);
3734 if ( (ir->eConstrAlg == econtLINCS) && bConstr)
3736 /* If we have Lincs constraints: */
3737 if (ir->eI == eiMD && ir->etc == etcNO &&
3738 ir->eConstrAlg == econtLINCS && ir->nLincsIter == 1)
3740 sprintf(warn_buf, "For energy conservation with LINCS, lincs_iter should be 2 or larger.\n");
3741 warning_note(wi, warn_buf);
3744 if ((ir->eI == eiCG || ir->eI == eiLBFGS) && (ir->nProjOrder < 8))
3746 sprintf(warn_buf, "For accurate %s with LINCS constraints, lincs-order should be 8 or more.", ei_names[ir->eI]);
3747 warning_note(wi, warn_buf);
3749 if (ir->epc == epcMTTK)
3751 warning_error(wi, "MTTK not compatible with lincs -- use shake instead.");
3755 if (ir->LincsWarnAngle > 90.0)
3757 sprintf(warn_buf, "lincs-warnangle can not be larger than 90 degrees, setting it to 90.\n");
3758 warning(wi, warn_buf);
3759 ir->LincsWarnAngle = 90.0;
3762 if (ir->ePBC != epbcNONE)
3764 if (ir->nstlist == 0)
3766 warning(wi, "With nstlist=0 atoms are only put into the box at step 0, therefore drifting atoms might cause the simulation to crash.");
3768 bTWIN = (ir->rlistlong > ir->rlist);
3769 if (ir->ns_type == ensGRID)
3771 if (sqr(ir->rlistlong) >= max_cutoff2(ir->ePBC, box))
3773 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",
3774 bTWIN ? (ir->rcoulomb == ir->rlistlong ? "rcoulomb" : "rvdw") : "rlist");
3775 warning_error(wi, warn_buf);
3780 min_size = min(box[XX][XX], min(box[YY][YY], box[ZZ][ZZ]));
3781 if (2*ir->rlistlong >= min_size)
3783 sprintf(warn_buf, "ERROR: One of the box lengths is smaller than twice the cut-off length. Increase the box size or decrease rlist.");
3784 warning_error(wi, warn_buf);
3787 fprintf(stderr, "Grid search might allow larger cut-off's than simple search with triclinic boxes.");
3794 void check_chargegroup_radii(const gmx_mtop_t *mtop, const t_inputrec *ir,
3798 real rvdw1, rvdw2, rcoul1, rcoul2;
3799 char warn_buf[STRLEN];
3801 calc_chargegroup_radii(mtop, x, &rvdw1, &rvdw2, &rcoul1, &rcoul2);
3805 printf("Largest charge group radii for Van der Waals: %5.3f, %5.3f nm\n",
3810 printf("Largest charge group radii for Coulomb: %5.3f, %5.3f nm\n",
3816 if (rvdw1 + rvdw2 > ir->rlist ||
3817 rcoul1 + rcoul2 > ir->rlist)
3820 "The sum of the two largest charge group radii (%f) "
3821 "is larger than rlist (%f)\n",
3822 max(rvdw1+rvdw2, rcoul1+rcoul2), ir->rlist);
3823 warning(wi, warn_buf);
3827 /* Here we do not use the zero at cut-off macro,
3828 * since user defined interactions might purposely
3829 * not be zero at the cut-off.
3831 if ((EVDW_IS_ZERO_AT_CUTOFF(ir->vdwtype) ||
3832 ir->vdw_modifier != eintmodNONE) &&
3833 rvdw1 + rvdw2 > ir->rlistlong - ir->rvdw)
3835 sprintf(warn_buf, "The sum of the two largest charge group "
3836 "radii (%f) is larger than %s (%f) - rvdw (%f).\n"
3837 "With exact cut-offs, better performance can be "
3838 "obtained with cutoff-scheme = %s, because it "
3839 "does not use charge groups at all.",
3841 ir->rlistlong > ir->rlist ? "rlistlong" : "rlist",
3842 ir->rlistlong, ir->rvdw,
3843 ecutscheme_names[ecutsVERLET]);
3846 warning(wi, warn_buf);
3850 warning_note(wi, warn_buf);
3853 if ((EEL_IS_ZERO_AT_CUTOFF(ir->coulombtype) ||
3854 ir->coulomb_modifier != eintmodNONE) &&
3855 rcoul1 + rcoul2 > ir->rlistlong - ir->rcoulomb)
3857 sprintf(warn_buf, "The sum of the two largest charge group radii (%f) is larger than %s (%f) - rcoulomb (%f).\n"
3858 "With exact cut-offs, better performance can be obtained with cutoff-scheme = %s, because it does not use charge groups at all.",
3860 ir->rlistlong > ir->rlist ? "rlistlong" : "rlist",
3861 ir->rlistlong, ir->rcoulomb,
3862 ecutscheme_names[ecutsVERLET]);
3865 warning(wi, warn_buf);
3869 warning_note(wi, warn_buf);