2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
49 #include "gmx_fatal.h"
62 #include "mtop_util.h"
63 #include "chargegroup.h"
68 #define MAXLAMBDAS 1024
70 /* Resource parameters
71 * Do not change any of these until you read the instruction
72 * in readinp.h. Some cpp's do not take spaces after the backslash
73 * (like the c-shell), which will give you a very weird compiler
77 static char tcgrps[STRLEN], tau_t[STRLEN], ref_t[STRLEN],
78 acc[STRLEN], accgrps[STRLEN], freeze[STRLEN], frdim[STRLEN],
79 energy[STRLEN], user1[STRLEN], user2[STRLEN], vcm[STRLEN], xtc_grps[STRLEN],
80 couple_moltype[STRLEN], orirefitgrp[STRLEN], egptable[STRLEN], egpexcl[STRLEN],
81 wall_atomtype[STRLEN], wall_density[STRLEN], deform[STRLEN], QMMM[STRLEN];
82 static char fep_lambda[efptNR][STRLEN];
83 static char lambda_weights[STRLEN];
84 static char **pull_grp;
85 static char **rot_grp;
86 static char anneal[STRLEN], anneal_npoints[STRLEN],
87 anneal_time[STRLEN], anneal_temp[STRLEN];
88 static char QMmethod[STRLEN], QMbasis[STRLEN], QMcharge[STRLEN], QMmult[STRLEN],
89 bSH[STRLEN], CASorbitals[STRLEN], CASelectrons[STRLEN], SAon[STRLEN],
90 SAoff[STRLEN], SAsteps[STRLEN], bTS[STRLEN], bOPT[STRLEN];
91 static char efield_x[STRLEN], efield_xt[STRLEN], efield_y[STRLEN],
92 efield_yt[STRLEN], efield_z[STRLEN], efield_zt[STRLEN];
95 egrptpALL, /* All particles have to be a member of a group. */
96 egrptpALL_GENREST, /* A rest group with name is generated for particles *
97 * that are not part of any group. */
98 egrptpPART, /* As egrptpALL_GENREST, but no name is generated *
99 * for the rest group. */
100 egrptpONE /* Merge all selected groups into one group, *
101 * make a rest group for the remaining particles. */
104 static const char *constraints[eshNR+1] = {
105 "none", "h-bonds", "all-bonds", "h-angles", "all-angles", NULL
108 static const char *couple_lam[ecouplamNR+1] = {
109 "vdw-q", "vdw", "q", "none", NULL
112 void init_ir(t_inputrec *ir, t_gromppopts *opts)
114 snew(opts->include, STRLEN);
115 snew(opts->define, STRLEN);
116 snew(ir->fepvals, 1);
117 snew(ir->expandedvals, 1);
118 snew(ir->simtempvals, 1);
121 static void GetSimTemps(int ntemps, t_simtemp *simtemp, double *temperature_lambdas)
126 for (i = 0; i < ntemps; i++)
128 /* simple linear scaling -- allows more control */
129 if (simtemp->eSimTempScale == esimtempLINEAR)
131 simtemp->temperatures[i] = simtemp->simtemp_low + (simtemp->simtemp_high-simtemp->simtemp_low)*temperature_lambdas[i];
133 else if (simtemp->eSimTempScale == esimtempGEOMETRIC) /* should give roughly equal acceptance for constant heat capacity . . . */
135 simtemp->temperatures[i] = simtemp->simtemp_low * pow(simtemp->simtemp_high/simtemp->simtemp_low, (1.0*i)/(ntemps-1));
137 else if (simtemp->eSimTempScale == esimtempEXPONENTIAL)
139 simtemp->temperatures[i] = simtemp->simtemp_low + (simtemp->simtemp_high-simtemp->simtemp_low)*((exp(temperature_lambdas[i])-1)/(exp(1.0)-1));
144 sprintf(errorstr, "eSimTempScale=%d not defined", simtemp->eSimTempScale);
145 gmx_fatal(FARGS, errorstr);
152 static void _low_check(gmx_bool b, char *s, warninp_t wi)
156 warning_error(wi, s);
160 static void check_nst(const char *desc_nst, int nst,
161 const char *desc_p, int *p,
166 if (*p > 0 && *p % nst != 0)
168 /* Round up to the next multiple of nst */
169 *p = ((*p)/nst + 1)*nst;
170 sprintf(buf, "%s should be a multiple of %s, changing %s to %d\n",
171 desc_p, desc_nst, desc_p, *p);
176 static gmx_bool ir_NVE(const t_inputrec *ir)
178 return ((ir->eI == eiMD || EI_VV(ir->eI)) && ir->etc == etcNO);
181 static int lcd(int n1, int n2)
186 for (i = 2; (i <= n1 && i <= n2); i++)
188 if (n1 % i == 0 && n2 % i == 0)
197 static void process_interaction_modifier(const t_inputrec *ir, int *eintmod)
199 if (*eintmod == eintmodPOTSHIFT_VERLET)
201 if (ir->cutoff_scheme == ecutsVERLET)
203 *eintmod = eintmodPOTSHIFT;
207 *eintmod = eintmodNONE;
212 void check_ir(const char *mdparin, t_inputrec *ir, t_gromppopts *opts,
214 /* Check internal consistency */
216 /* Strange macro: first one fills the err_buf, and then one can check
217 * the condition, which will print the message and increase the error
220 #define CHECK(b) _low_check(b, err_buf, wi)
221 char err_buf[256], warn_buf[STRLEN];
227 t_lambda *fep = ir->fepvals;
228 t_expanded *expand = ir->expandedvals;
230 set_warning_line(wi, mdparin, -1);
232 /* BASIC CUT-OFF STUFF */
233 if (ir->rcoulomb < 0)
235 warning_error(wi, "rcoulomb should be >= 0");
239 warning_error(wi, "rvdw should be >= 0");
242 !(ir->cutoff_scheme == ecutsVERLET && ir->verletbuf_tol > 0))
244 warning_error(wi, "rlist should be >= 0");
247 process_interaction_modifier(ir, &ir->coulomb_modifier);
248 process_interaction_modifier(ir, &ir->vdw_modifier);
250 if (ir->cutoff_scheme == ecutsGROUP)
252 /* BASIC CUT-OFF STUFF */
253 if (ir->rlist == 0 ||
254 !((EEL_MIGHT_BE_ZERO_AT_CUTOFF(ir->coulombtype) && ir->rcoulomb > ir->rlist) ||
255 (EVDW_MIGHT_BE_ZERO_AT_CUTOFF(ir->vdwtype) && ir->rvdw > ir->rlist)))
257 /* No switched potential and/or no twin-range:
258 * we can set the long-range cut-off to the maximum of the other cut-offs.
260 ir->rlistlong = max_cutoff(ir->rlist, max_cutoff(ir->rvdw, ir->rcoulomb));
262 else if (ir->rlistlong < 0)
264 ir->rlistlong = max_cutoff(ir->rlist, max_cutoff(ir->rvdw, ir->rcoulomb));
265 sprintf(warn_buf, "rlistlong was not set, setting it to %g (no buffer)",
267 warning(wi, warn_buf);
269 if (ir->rlistlong == 0 && ir->ePBC != epbcNONE)
271 warning_error(wi, "Can not have an infinite cut-off with PBC");
273 if (ir->rlistlong > 0 && (ir->rlist == 0 || ir->rlistlong < ir->rlist))
275 warning_error(wi, "rlistlong can not be shorter than rlist");
277 if (IR_TWINRANGE(*ir) && ir->nstlist <= 0)
279 warning_error(wi, "Can not have nstlist<=0 with twin-range interactions");
283 if (ir->rlistlong == ir->rlist)
287 else if (ir->rlistlong > ir->rlist && ir->nstcalclr == 0)
289 warning_error(wi, "With different cutoffs for electrostatics and VdW, nstcalclr must be -1 or a positive number");
292 if (ir->cutoff_scheme == ecutsVERLET)
296 /* Normal Verlet type neighbor-list, currently only limited feature support */
297 if (inputrec2nboundeddim(ir) < 3)
299 warning_error(wi, "With Verlet lists only full pbc or pbc=xy with walls is supported");
301 if (ir->rcoulomb != ir->rvdw)
303 warning_error(wi, "With Verlet lists rcoulomb!=rvdw is not supported");
305 if (ir->vdwtype != evdwCUT)
307 warning_error(wi, "With Verlet lists only cut-off LJ interactions are supported");
309 if (!(ir->coulombtype == eelCUT ||
310 (EEL_RF(ir->coulombtype) && ir->coulombtype != eelRF_NEC) ||
311 EEL_PME(ir->coulombtype) || ir->coulombtype == eelEWALD))
313 warning_error(wi, "With Verlet lists only cut-off, reaction-field, PME and Ewald electrostatics are supported");
316 if (ir->nstlist <= 0)
318 warning_error(wi, "With Verlet lists nstlist should be larger than 0");
321 if (ir->nstlist < 10)
323 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.");
326 rc_max = max(ir->rvdw, ir->rcoulomb);
328 if (ir->verletbuf_tol <= 0)
330 if (ir->verletbuf_tol == 0)
332 warning_error(wi, "Can not have Verlet buffer tolerance of exactly 0");
335 if (ir->rlist < rc_max)
337 warning_error(wi, "With verlet lists rlist can not be smaller than rvdw or rcoulomb");
340 if (ir->rlist == rc_max && ir->nstlist > 1)
342 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.");
347 if (ir->rlist > rc_max)
349 warning_note(wi, "You have set rlist larger than the interaction cut-off, but you also have verlet-buffer-tolerance > 0. Will set rlist using verlet-buffer-tolerance.");
352 if (ir->nstlist == 1)
354 /* No buffer required */
359 if (EI_DYNAMICS(ir->eI))
361 if (EI_MD(ir->eI) && ir->etc == etcNO)
363 warning_error(wi, "Temperature coupling is required for calculating rlist using the energy tolerance with verlet-buffer-tolerance > 0. Either use temperature coupling or set rlist yourself together with verlet-buffer-tolerance = -1.");
366 if (inputrec2nboundeddim(ir) < 3)
368 warning_error(wi, "The box volume is required for calculating rlist from the energy drift with verlet-buffer-tolerance > 0. You are using at least one unbounded dimension, so no volume can be computed. Either use a finite box, or set rlist yourself together with verlet-buffer-tolerance = -1.");
370 /* Set rlist temporarily so we can continue processing */
375 /* Set the buffer to 5% of the cut-off */
376 ir->rlist = 1.05*rc_max;
381 /* No twin-range calculations with Verlet lists */
382 ir->rlistlong = ir->rlist;
385 if (ir->nstcalclr == -1)
387 /* if rlist=rlistlong, this will later be changed to nstcalclr=0 */
388 ir->nstcalclr = ir->nstlist;
390 else if (ir->nstcalclr > 0)
392 if (ir->nstlist > 0 && (ir->nstlist % ir->nstcalclr != 0))
394 warning_error(wi, "nstlist must be evenly divisible by nstcalclr. Use nstcalclr = -1 to automatically follow nstlist");
397 else if (ir->nstcalclr < -1)
399 warning_error(wi, "nstcalclr must be a positive number (divisor of nstcalclr), or -1 to follow nstlist.");
402 if (EEL_PME(ir->coulombtype) && ir->rcoulomb > ir->rvdw && ir->nstcalclr > 1)
404 warning_error(wi, "When used with PME, the long-range component of twin-range interactions must be updated every step (nstcalclr)");
407 /* GENERAL INTEGRATOR STUFF */
408 if (!(ir->eI == eiMD || EI_VV(ir->eI)))
412 if (ir->eI == eiVVAK)
414 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]);
415 warning_note(wi, warn_buf);
417 if (!EI_DYNAMICS(ir->eI))
421 if (EI_DYNAMICS(ir->eI))
423 if (ir->nstcalcenergy < 0)
425 ir->nstcalcenergy = ir_optimal_nstcalcenergy(ir);
426 if (ir->nstenergy != 0 && ir->nstenergy < ir->nstcalcenergy)
428 /* nstcalcenergy larger than nstener does not make sense.
429 * We ideally want nstcalcenergy=nstener.
433 ir->nstcalcenergy = lcd(ir->nstenergy, ir->nstlist);
437 ir->nstcalcenergy = ir->nstenergy;
441 else if ( (ir->nstenergy > 0 && ir->nstcalcenergy > ir->nstenergy) ||
442 (ir->efep != efepNO && ir->fepvals->nstdhdl > 0 &&
443 (ir->nstcalcenergy > ir->fepvals->nstdhdl) ) )
446 const char *nsten = "nstenergy";
447 const char *nstdh = "nstdhdl";
448 const char *min_name = nsten;
449 int min_nst = ir->nstenergy;
451 /* find the smallest of ( nstenergy, nstdhdl ) */
452 if (ir->efep != efepNO && ir->fepvals->nstdhdl > 0 &&
453 (ir->fepvals->nstdhdl < ir->nstenergy) )
455 min_nst = ir->fepvals->nstdhdl;
458 /* If the user sets nstenergy small, we should respect that */
460 "Setting nstcalcenergy (%d) equal to %s (%d)",
461 ir->nstcalcenergy, min_name, min_nst);
462 warning_note(wi, warn_buf);
463 ir->nstcalcenergy = min_nst;
466 if (ir->epc != epcNO)
468 if (ir->nstpcouple < 0)
470 ir->nstpcouple = ir_optimal_nstpcouple(ir);
473 if (IR_TWINRANGE(*ir))
475 check_nst("nstlist", ir->nstlist,
476 "nstcalcenergy", &ir->nstcalcenergy, wi);
477 if (ir->epc != epcNO)
479 check_nst("nstlist", ir->nstlist,
480 "nstpcouple", &ir->nstpcouple, wi);
484 if (ir->nstcalcenergy > 0)
486 if (ir->efep != efepNO)
488 /* nstdhdl should be a multiple of nstcalcenergy */
489 check_nst("nstcalcenergy", ir->nstcalcenergy,
490 "nstdhdl", &ir->fepvals->nstdhdl, wi);
491 /* nstexpanded should be a multiple of nstcalcenergy */
492 check_nst("nstcalcenergy", ir->nstcalcenergy,
493 "nstexpanded", &ir->expandedvals->nstexpanded, wi);
495 /* for storing exact averages nstenergy should be
496 * a multiple of nstcalcenergy
498 check_nst("nstcalcenergy", ir->nstcalcenergy,
499 "nstenergy", &ir->nstenergy, wi);
504 if ((EI_SD(ir->eI) || ir->eI == eiBD) &&
505 ir->bContinuation && ir->ld_seed != -1)
507 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)");
513 sprintf(err_buf, "TPI only works with pbc = %s", epbc_names[epbcXYZ]);
514 CHECK(ir->ePBC != epbcXYZ);
515 sprintf(err_buf, "TPI only works with ns = %s", ens_names[ensGRID]);
516 CHECK(ir->ns_type != ensGRID);
517 sprintf(err_buf, "with TPI nstlist should be larger than zero");
518 CHECK(ir->nstlist <= 0);
519 sprintf(err_buf, "TPI does not work with full electrostatics other than PME");
520 CHECK(EEL_FULL(ir->coulombtype) && !EEL_PME(ir->coulombtype));
524 if ( (opts->nshake > 0) && (opts->bMorse) )
527 "Using morse bond-potentials while constraining bonds is useless");
528 warning(wi, warn_buf);
531 if ((EI_SD(ir->eI) || ir->eI == eiBD) &&
532 ir->bContinuation && ir->ld_seed != -1)
534 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)");
536 /* verify simulated tempering options */
540 gmx_bool bAllTempZero = TRUE;
541 for (i = 0; i < fep->n_lambda; i++)
543 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]);
544 CHECK((fep->all_lambda[efptTEMPERATURE][i] < 0) || (fep->all_lambda[efptTEMPERATURE][i] > 1));
545 if (fep->all_lambda[efptTEMPERATURE][i] > 0)
547 bAllTempZero = FALSE;
550 sprintf(err_buf, "if simulated tempering is on, temperature-lambdas may not be all zero");
551 CHECK(bAllTempZero == TRUE);
553 sprintf(err_buf, "Simulated tempering is currently only compatible with md-vv");
554 CHECK(ir->eI != eiVV);
556 /* check compatability of the temperature coupling with simulated tempering */
558 if (ir->etc == etcNOSEHOOVER)
560 sprintf(warn_buf, "Nose-Hoover based temperature control such as [%s] my not be entirelyconsistent with simulated tempering", etcoupl_names[ir->etc]);
561 warning_note(wi, warn_buf);
564 /* check that the temperatures make sense */
566 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);
567 CHECK(ir->simtempvals->simtemp_high <= ir->simtempvals->simtemp_low);
569 sprintf(err_buf, "Higher simulated tempering temperature (%g) must be >= zero", ir->simtempvals->simtemp_high);
570 CHECK(ir->simtempvals->simtemp_high <= 0);
572 sprintf(err_buf, "Lower simulated tempering temperature (%g) must be >= zero", ir->simtempvals->simtemp_low);
573 CHECK(ir->simtempvals->simtemp_low <= 0);
576 /* verify free energy options */
578 if (ir->efep != efepNO)
581 sprintf(err_buf, "The soft-core power is %d and can only be 1 or 2",
583 CHECK(fep->sc_alpha != 0 && fep->sc_power != 1 && fep->sc_power != 2);
585 sprintf(err_buf, "The soft-core sc-r-power is %d and can only be 6 or 48",
586 (int)fep->sc_r_power);
587 CHECK(fep->sc_alpha != 0 && fep->sc_r_power != 6.0 && fep->sc_r_power != 48.0);
589 sprintf(err_buf, "Can't use postive delta-lambda (%g) if initial state/lambda does not start at zero", fep->delta_lambda);
590 CHECK(fep->delta_lambda > 0 && ((fep->init_fep_state > 0) || (fep->init_lambda > 0)));
592 sprintf(err_buf, "Can't use postive delta-lambda (%g) with expanded ensemble simulations", fep->delta_lambda);
593 CHECK(fep->delta_lambda > 0 && (ir->efep == efepEXPANDED));
595 sprintf(err_buf, "Can only use expanded ensemble with md-vv for now; should be supported for other integrators in 5.0");
596 CHECK(!(EI_VV(ir->eI)) && (ir->efep == efepEXPANDED));
598 sprintf(err_buf, "Free-energy not implemented for Ewald");
599 CHECK(ir->coulombtype == eelEWALD);
601 /* check validty of lambda inputs */
602 if (fep->n_lambda == 0)
604 /* Clear output in case of no states:*/
605 sprintf(err_buf, "init-lambda-state set to %d: no lambda states are defined.", fep->init_fep_state);
606 CHECK((fep->init_fep_state >= 0) && (fep->n_lambda == 0));
610 sprintf(err_buf, "initial thermodynamic state %d does not exist, only goes to %d", fep->init_fep_state, fep->n_lambda-1);
611 CHECK((fep->init_fep_state >= fep->n_lambda));
614 sprintf(err_buf, "Lambda state must be set, either with init-lambda-state or with init-lambda");
615 CHECK((fep->init_fep_state < 0) && (fep->init_lambda < 0));
617 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",
618 fep->init_lambda, fep->init_fep_state);
619 CHECK((fep->init_fep_state >= 0) && (fep->init_lambda >= 0));
623 if ((fep->init_lambda >= 0) && (fep->delta_lambda == 0))
627 for (i = 0; i < efptNR; i++)
629 if (fep->separate_dvdl[i])
634 if (n_lambda_terms > 1)
636 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.");
637 warning(wi, warn_buf);
640 if (n_lambda_terms < 2 && fep->n_lambda > 0)
643 "init-lambda is deprecated for setting lambda state (except for slow growth). Use init-lambda-state instead.");
647 for (j = 0; j < efptNR; j++)
649 for (i = 0; i < fep->n_lambda; i++)
651 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]);
652 CHECK((fep->all_lambda[j][i] < 0) || (fep->all_lambda[j][i] > 1));
656 if ((fep->sc_alpha > 0) && (!fep->bScCoul))
658 for (i = 0; i < fep->n_lambda; i++)
660 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],
661 fep->all_lambda[efptCOUL][i]);
662 CHECK((fep->sc_alpha > 0) &&
663 (((fep->all_lambda[efptCOUL][i] > 0.0) &&
664 (fep->all_lambda[efptCOUL][i] < 1.0)) &&
665 ((fep->all_lambda[efptVDW][i] > 0.0) &&
666 (fep->all_lambda[efptVDW][i] < 1.0))));
670 if ((fep->bScCoul) && (EEL_PME(ir->coulombtype)))
672 real sigma, lambda, r_sc;
675 /* Maximum estimate for A and B charges equal with lambda power 1 */
677 r_sc = pow(lambda*fep->sc_alpha*pow(sigma/ir->rcoulomb, fep->sc_r_power) + 1.0, 1.0/fep->sc_r_power);
678 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.",
680 sigma, lambda, r_sc - 1.0, ir->ewald_rtol);
681 warning_note(wi, warn_buf);
684 /* Free Energy Checks -- In an ideal world, slow growth and FEP would
685 be treated differently, but that's the next step */
687 for (i = 0; i < efptNR; i++)
689 for (j = 0; j < fep->n_lambda; j++)
691 sprintf(err_buf, "%s[%d] must be between 0 and 1", efpt_names[i], j);
692 CHECK((fep->all_lambda[i][j] < 0) || (fep->all_lambda[i][j] > 1));
697 if ((ir->bSimTemp) || (ir->efep == efepEXPANDED))
700 expand = ir->expandedvals;
702 /* checking equilibration of weights inputs for validity */
704 sprintf(err_buf, "weight-equil-number-all-lambda (%d) is ignored if lmc-weights-equil is not equal to %s",
705 expand->equil_n_at_lam, elmceq_names[elmceqNUMATLAM]);
706 CHECK((expand->equil_n_at_lam > 0) && (expand->elmceq != elmceqNUMATLAM));
708 sprintf(err_buf, "weight-equil-number-samples (%d) is ignored if lmc-weights-equil is not equal to %s",
709 expand->equil_samples, elmceq_names[elmceqSAMPLES]);
710 CHECK((expand->equil_samples > 0) && (expand->elmceq != elmceqSAMPLES));
712 sprintf(err_buf, "weight-equil-number-steps (%d) is ignored if lmc-weights-equil is not equal to %s",
713 expand->equil_steps, elmceq_names[elmceqSTEPS]);
714 CHECK((expand->equil_steps > 0) && (expand->elmceq != elmceqSTEPS));
716 sprintf(err_buf, "weight-equil-wl-delta (%d) is ignored if lmc-weights-equil is not equal to %s",
717 expand->equil_samples, elmceq_names[elmceqWLDELTA]);
718 CHECK((expand->equil_wl_delta > 0) && (expand->elmceq != elmceqWLDELTA));
720 sprintf(err_buf, "weight-equil-count-ratio (%f) is ignored if lmc-weights-equil is not equal to %s",
721 expand->equil_ratio, elmceq_names[elmceqRATIO]);
722 CHECK((expand->equil_ratio > 0) && (expand->elmceq != elmceqRATIO));
724 sprintf(err_buf, "weight-equil-number-all-lambda (%d) must be a positive integer if lmc-weights-equil=%s",
725 expand->equil_n_at_lam, elmceq_names[elmceqNUMATLAM]);
726 CHECK((expand->equil_n_at_lam <= 0) && (expand->elmceq == elmceqNUMATLAM));
728 sprintf(err_buf, "weight-equil-number-samples (%d) must be a positive integer if lmc-weights-equil=%s",
729 expand->equil_samples, elmceq_names[elmceqSAMPLES]);
730 CHECK((expand->equil_samples <= 0) && (expand->elmceq == elmceqSAMPLES));
732 sprintf(err_buf, "weight-equil-number-steps (%d) must be a positive integer if lmc-weights-equil=%s",
733 expand->equil_steps, elmceq_names[elmceqSTEPS]);
734 CHECK((expand->equil_steps <= 0) && (expand->elmceq == elmceqSTEPS));
736 sprintf(err_buf, "weight-equil-wl-delta (%f) must be > 0 if lmc-weights-equil=%s",
737 expand->equil_wl_delta, elmceq_names[elmceqWLDELTA]);
738 CHECK((expand->equil_wl_delta <= 0) && (expand->elmceq == elmceqWLDELTA));
740 sprintf(err_buf, "weight-equil-count-ratio (%f) must be > 0 if lmc-weights-equil=%s",
741 expand->equil_ratio, elmceq_names[elmceqRATIO]);
742 CHECK((expand->equil_ratio <= 0) && (expand->elmceq == elmceqRATIO));
744 sprintf(err_buf, "lmc-weights-equil=%s only possible when lmc-stats = %s or lmc-stats %s",
745 elmceq_names[elmceqWLDELTA], elamstats_names[elamstatsWL], elamstats_names[elamstatsWWL]);
746 CHECK((expand->elmceq == elmceqWLDELTA) && (!EWL(expand->elamstats)));
748 sprintf(err_buf, "lmc-repeats (%d) must be greater than 0", expand->lmc_repeats);
749 CHECK((expand->lmc_repeats <= 0));
750 sprintf(err_buf, "minimum-var-min (%d) must be greater than 0", expand->minvarmin);
751 CHECK((expand->minvarmin <= 0));
752 sprintf(err_buf, "weight-c-range (%d) must be greater or equal to 0", expand->c_range);
753 CHECK((expand->c_range < 0));
754 sprintf(err_buf, "init-lambda-state (%d) must be zero if lmc-forced-nstart (%d)> 0 and lmc-move != 'no'",
755 fep->init_fep_state, expand->lmc_forced_nstart);
756 CHECK((fep->init_fep_state != 0) && (expand->lmc_forced_nstart > 0) && (expand->elmcmove != elmcmoveNO));
757 sprintf(err_buf, "lmc-forced-nstart (%d) must not be negative", expand->lmc_forced_nstart);
758 CHECK((expand->lmc_forced_nstart < 0));
759 sprintf(err_buf, "init-lambda-state (%d) must be in the interval [0,number of lambdas)", fep->init_fep_state);
760 CHECK((fep->init_fep_state < 0) || (fep->init_fep_state >= fep->n_lambda));
762 sprintf(err_buf, "init-wl-delta (%f) must be greater than or equal to 0", expand->init_wl_delta);
763 CHECK((expand->init_wl_delta < 0));
764 sprintf(err_buf, "wl-ratio (%f) must be between 0 and 1", expand->wl_ratio);
765 CHECK((expand->wl_ratio <= 0) || (expand->wl_ratio >= 1));
766 sprintf(err_buf, "wl-scale (%f) must be between 0 and 1", expand->wl_scale);
767 CHECK((expand->wl_scale <= 0) || (expand->wl_scale >= 1));
769 /* if there is no temperature control, we need to specify an MC temperature */
770 sprintf(err_buf, "If there is no temperature control, and lmc-mcmove!= 'no',mc_temperature must be set to a positive number");
771 if (expand->nstTij > 0)
773 sprintf(err_buf, "nst-transition-matrix (%d) must be an integer multiple of nstlog (%d)",
774 expand->nstTij, ir->nstlog);
775 CHECK((mod(expand->nstTij, ir->nstlog) != 0));
780 sprintf(err_buf, "walls only work with pbc=%s", epbc_names[epbcXY]);
781 CHECK(ir->nwall && ir->ePBC != epbcXY);
784 if (ir->ePBC != epbcXYZ && ir->nwall != 2)
786 if (ir->ePBC == epbcNONE)
788 if (ir->epc != epcNO)
790 warning(wi, "Turning off pressure coupling for vacuum system");
796 sprintf(err_buf, "Can not have pressure coupling with pbc=%s",
797 epbc_names[ir->ePBC]);
798 CHECK(ir->epc != epcNO);
800 sprintf(err_buf, "Can not have Ewald with pbc=%s", epbc_names[ir->ePBC]);
801 CHECK(EEL_FULL(ir->coulombtype));
803 sprintf(err_buf, "Can not have dispersion correction with pbc=%s",
804 epbc_names[ir->ePBC]);
805 CHECK(ir->eDispCorr != edispcNO);
808 if (ir->rlist == 0.0)
810 sprintf(err_buf, "can only have neighborlist cut-off zero (=infinite)\n"
811 "with coulombtype = %s or coulombtype = %s\n"
812 "without periodic boundary conditions (pbc = %s) and\n"
813 "rcoulomb and rvdw set to zero",
814 eel_names[eelCUT], eel_names[eelUSER], epbc_names[epbcNONE]);
815 CHECK(((ir->coulombtype != eelCUT) && (ir->coulombtype != eelUSER)) ||
816 (ir->ePBC != epbcNONE) ||
817 (ir->rcoulomb != 0.0) || (ir->rvdw != 0.0));
821 warning_error(wi, "Can not have heuristic neighborlist updates without cut-off");
825 warning_note(wi, "Simulating without cut-offs is usually (slightly) faster with nstlist=0, nstype=simple and particle decomposition");
830 if (ir->nstcomm == 0)
832 ir->comm_mode = ecmNO;
834 if (ir->comm_mode != ecmNO)
838 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");
839 ir->nstcomm = abs(ir->nstcomm);
842 if (ir->nstcalcenergy > 0 && ir->nstcomm < ir->nstcalcenergy)
844 warning_note(wi, "nstcomm < nstcalcenergy defeats the purpose of nstcalcenergy, setting nstcomm to nstcalcenergy");
845 ir->nstcomm = ir->nstcalcenergy;
848 if (ir->comm_mode == ecmANGULAR)
850 sprintf(err_buf, "Can not remove the rotation around the center of mass with periodic molecules");
851 CHECK(ir->bPeriodicMols);
852 if (ir->ePBC != epbcNONE)
854 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).");
859 if (EI_STATE_VELOCITY(ir->eI) && ir->ePBC == epbcNONE && ir->comm_mode != ecmANGULAR)
861 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.");
864 sprintf(err_buf, "Twin-range neighbour searching (NS) with simple NS"
865 " algorithm not implemented");
866 CHECK(((ir->rcoulomb > ir->rlist) || (ir->rvdw > ir->rlist))
867 && (ir->ns_type == ensSIMPLE));
869 /* TEMPERATURE COUPLING */
870 if (ir->etc == etcYES)
872 ir->etc = etcBERENDSEN;
873 warning_note(wi, "Old option for temperature coupling given: "
874 "changing \"yes\" to \"Berendsen\"\n");
877 if ((ir->etc == etcNOSEHOOVER) || (ir->epc == epcMTTK))
879 if (ir->opts.nhchainlength < 1)
881 sprintf(warn_buf, "number of Nose-Hoover chains (currently %d) cannot be less than 1,reset to 1\n", ir->opts.nhchainlength);
882 ir->opts.nhchainlength = 1;
883 warning(wi, warn_buf);
886 if (ir->etc == etcNOSEHOOVER && !EI_VV(ir->eI) && ir->opts.nhchainlength > 1)
888 warning_note(wi, "leapfrog does not yet support Nose-Hoover chains, nhchainlength reset to 1");
889 ir->opts.nhchainlength = 1;
894 ir->opts.nhchainlength = 0;
897 if (ir->eI == eiVVAK)
899 sprintf(err_buf, "%s implemented primarily for validation, and requires nsttcouple = 1 and nstpcouple = 1.",
901 CHECK((ir->nsttcouple != 1) || (ir->nstpcouple != 1));
904 if (ETC_ANDERSEN(ir->etc))
906 sprintf(err_buf, "%s temperature control not supported for integrator %s.", etcoupl_names[ir->etc], ei_names[ir->eI]);
907 CHECK(!(EI_VV(ir->eI)));
909 for (i = 0; i < ir->opts.ngtc; i++)
911 sprintf(err_buf, "all tau_t must currently be equal using Andersen temperature control, violated for group %d", i);
912 CHECK(ir->opts.tau_t[0] != ir->opts.tau_t[i]);
913 sprintf(err_buf, "all tau_t must be postive using Andersen temperature control, tau_t[%d]=%10.6f",
914 i, ir->opts.tau_t[i]);
915 CHECK(ir->opts.tau_t[i] < 0);
917 if (ir->nstcomm > 0 && (ir->etc == etcANDERSEN))
919 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]);
920 warning_note(wi, warn_buf);
923 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]);
924 CHECK(ir->nstcomm > 1 && (ir->etc == etcANDERSEN));
926 for (i = 0; i < ir->opts.ngtc; i++)
928 int nsteps = (int)(ir->opts.tau_t[i]/ir->delta_t);
929 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);
930 CHECK((nsteps % ir->nstcomm) && (ir->etc == etcANDERSENMASSIVE));
933 if (ir->etc == etcBERENDSEN)
935 sprintf(warn_buf, "The %s thermostat does not generate the correct kinetic energy distribution. You might want to consider using the %s thermostat.",
936 ETCOUPLTYPE(ir->etc), ETCOUPLTYPE(etcVRESCALE));
937 warning_note(wi, warn_buf);
940 if ((ir->etc == etcNOSEHOOVER || ETC_ANDERSEN(ir->etc))
941 && ir->epc == epcBERENDSEN)
943 sprintf(warn_buf, "Using Berendsen pressure coupling invalidates the "
944 "true ensemble for the thermostat");
945 warning(wi, warn_buf);
948 /* PRESSURE COUPLING */
949 if (ir->epc == epcISOTROPIC)
951 ir->epc = epcBERENDSEN;
952 warning_note(wi, "Old option for pressure coupling given: "
953 "changing \"Isotropic\" to \"Berendsen\"\n");
956 if (ir->epc != epcNO)
958 dt_pcoupl = ir->nstpcouple*ir->delta_t;
960 sprintf(err_buf, "tau-p must be > 0 instead of %g\n", ir->tau_p);
961 CHECK(ir->tau_p <= 0);
963 if (ir->tau_p/dt_pcoupl < pcouple_min_integration_steps(ir->epc))
965 sprintf(warn_buf, "For proper integration of the %s barostat, tau-p (%g) should be at least %d times larger than nstpcouple*dt (%g)",
966 EPCOUPLTYPE(ir->epc), ir->tau_p, pcouple_min_integration_steps(ir->epc), dt_pcoupl);
967 warning(wi, warn_buf);
970 sprintf(err_buf, "compressibility must be > 0 when using pressure"
971 " coupling %s\n", EPCOUPLTYPE(ir->epc));
972 CHECK(ir->compress[XX][XX] < 0 || ir->compress[YY][YY] < 0 ||
973 ir->compress[ZZ][ZZ] < 0 ||
974 (trace(ir->compress) == 0 && ir->compress[YY][XX] <= 0 &&
975 ir->compress[ZZ][XX] <= 0 && ir->compress[ZZ][YY] <= 0));
977 if (epcPARRINELLORAHMAN == ir->epc && opts->bGenVel)
980 "You are generating velocities so I am assuming you "
981 "are equilibrating a system. You are using "
982 "%s pressure coupling, but this can be "
983 "unstable for equilibration. If your system crashes, try "
984 "equilibrating first with Berendsen pressure coupling. If "
985 "you are not equilibrating the system, you can probably "
986 "ignore this warning.",
987 epcoupl_names[ir->epc]);
988 warning(wi, warn_buf);
996 if ((ir->epc != epcBERENDSEN) && (ir->epc != epcMTTK))
998 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.");
1003 /* ELECTROSTATICS */
1004 /* More checks are in triple check (grompp.c) */
1006 if (ir->coulombtype == eelSWITCH)
1008 sprintf(warn_buf, "coulombtype = %s is only for testing purposes and can lead to serious "
1009 "artifacts, advice: use coulombtype = %s",
1010 eel_names[ir->coulombtype],
1011 eel_names[eelRF_ZERO]);
1012 warning(wi, warn_buf);
1015 if (ir->epsilon_r != 1 && ir->implicit_solvent == eisGBSA)
1017 sprintf(warn_buf, "epsilon-r = %g with GB implicit solvent, will use this value for inner dielectric", ir->epsilon_r);
1018 warning_note(wi, warn_buf);
1021 if (EEL_RF(ir->coulombtype) && ir->epsilon_rf == 1 && ir->epsilon_r != 1)
1023 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);
1024 warning(wi, warn_buf);
1025 ir->epsilon_rf = ir->epsilon_r;
1026 ir->epsilon_r = 1.0;
1029 if (getenv("GALACTIC_DYNAMICS") == NULL)
1031 sprintf(err_buf, "epsilon-r must be >= 0 instead of %g\n", ir->epsilon_r);
1032 CHECK(ir->epsilon_r < 0);
1035 if (EEL_RF(ir->coulombtype))
1037 /* reaction field (at the cut-off) */
1039 if (ir->coulombtype == eelRF_ZERO)
1041 sprintf(warn_buf, "With coulombtype = %s, epsilon-rf must be 0, assuming you meant epsilon_rf=0",
1042 eel_names[ir->coulombtype]);
1043 CHECK(ir->epsilon_rf != 0);
1044 ir->epsilon_rf = 0.0;
1047 sprintf(err_buf, "epsilon-rf must be >= epsilon-r");
1048 CHECK((ir->epsilon_rf < ir->epsilon_r && ir->epsilon_rf != 0) ||
1049 (ir->epsilon_r == 0));
1050 if (ir->epsilon_rf == ir->epsilon_r)
1052 sprintf(warn_buf, "Using epsilon-rf = epsilon-r with %s does not make sense",
1053 eel_names[ir->coulombtype]);
1054 warning(wi, warn_buf);
1057 /* Allow rlist>rcoulomb for tabulated long range stuff. This just
1058 * means the interaction is zero outside rcoulomb, but it helps to
1059 * provide accurate energy conservation.
1061 if (EEL_MIGHT_BE_ZERO_AT_CUTOFF(ir->coulombtype))
1063 if (EEL_SWITCHED(ir->coulombtype))
1066 "With coulombtype = %s rcoulomb_switch must be < rcoulomb. Or, better: Use the potential modifier options!",
1067 eel_names[ir->coulombtype]);
1068 CHECK(ir->rcoulomb_switch >= ir->rcoulomb);
1071 else if (ir->coulombtype == eelCUT || EEL_RF(ir->coulombtype))
1073 if (ir->cutoff_scheme == ecutsGROUP && ir->coulomb_modifier == eintmodNONE)
1075 sprintf(err_buf, "With coulombtype = %s, rcoulomb should be >= rlist unless you use a potential modifier",
1076 eel_names[ir->coulombtype]);
1077 CHECK(ir->rlist > ir->rcoulomb);
1081 if (ir->coulombtype == eelSWITCH || ir->coulombtype == eelSHIFT ||
1082 ir->vdwtype == evdwSWITCH || ir->vdwtype == evdwSHIFT)
1085 "The switch/shift interaction settings are just for compatibility; you will get better "
1086 "performance from applying potential modifiers to your interactions!\n");
1087 warning_note(wi, warn_buf);
1090 if (ir->coulombtype == eelPMESWITCH)
1092 if (ir->rcoulomb_switch/ir->rcoulomb < 0.9499)
1094 sprintf(warn_buf, "The switching range for %s should be 5%% or less, energy conservation will be good anyhow, since ewald_rtol = %g",
1095 eel_names[ir->coulombtype],
1097 warning(wi, warn_buf);
1101 if (EEL_FULL(ir->coulombtype))
1103 if (ir->coulombtype == eelPMESWITCH || ir->coulombtype == eelPMEUSER ||
1104 ir->coulombtype == eelPMEUSERSWITCH)
1106 sprintf(err_buf, "With coulombtype = %s, rcoulomb must be <= rlist",
1107 eel_names[ir->coulombtype]);
1108 CHECK(ir->rcoulomb > ir->rlist);
1110 else if (ir->cutoff_scheme == ecutsGROUP && ir->coulomb_modifier == eintmodNONE)
1112 if (ir->coulombtype == eelPME || ir->coulombtype == eelP3M_AD)
1115 "With coulombtype = %s (without modifier), rcoulomb must be equal to rlist,\n"
1116 "or rlistlong if nstcalclr=1. For optimal energy conservation,consider using\n"
1117 "a potential modifier.", eel_names[ir->coulombtype]);
1118 if (ir->nstcalclr == 1)
1120 CHECK(ir->rcoulomb != ir->rlist && ir->rcoulomb != ir->rlistlong);
1124 CHECK(ir->rcoulomb != ir->rlist);
1130 if (EVDW_PME(ir->vdwtype))
1132 if (EVDW_MIGHT_BE_ZERO_AT_CUTOFF(ir->vdwtype))
1134 sprintf(err_buf, "With vdwtype = %s, rvdw must be <= rlist",
1135 evdw_names[ir->vdwtype]);
1136 CHECK(ir->rvdw > ir->rlist);
1141 "With vdwtype = %s, rvdw must be equal to rlist\n",
1142 evdw_names[ir->vdwtype]);
1143 CHECK(ir->rvdw != ir->rlist);
1147 if (EEL_PME(ir->coulombtype) || EVDW_PME(ir->vdwtype))
1149 if (ir->pme_order < 3)
1151 warning_error(wi, "pme-order can not be smaller than 3");
1155 if (ir->nwall == 2 && EEL_FULL(ir->coulombtype))
1157 if (ir->ewald_geometry == eewg3D)
1159 sprintf(warn_buf, "With pbc=%s you should use ewald-geometry=%s",
1160 epbc_names[ir->ePBC], eewg_names[eewg3DC]);
1161 warning(wi, warn_buf);
1163 /* This check avoids extra pbc coding for exclusion corrections */
1164 sprintf(err_buf, "wall-ewald-zfac should be >= 2");
1165 CHECK(ir->wall_ewald_zfac < 2);
1168 if (EVDW_SWITCHED(ir->vdwtype))
1170 sprintf(err_buf, "With vdwtype = %s rvdw-switch must be < rvdw. Or, better - use a potential modifier.",
1171 evdw_names[ir->vdwtype]);
1172 CHECK(ir->rvdw_switch >= ir->rvdw);
1174 else if (ir->vdwtype == evdwCUT)
1176 if (ir->cutoff_scheme == ecutsGROUP && ir->vdw_modifier == eintmodNONE)
1178 sprintf(err_buf, "With vdwtype = %s, rvdw must be >= rlist unless you use a potential modifier", evdw_names[ir->vdwtype]);
1179 CHECK(ir->rlist > ir->rvdw);
1182 if (ir->cutoff_scheme == ecutsGROUP)
1184 if (((ir->coulomb_modifier != eintmodNONE && ir->rcoulomb == ir->rlist) ||
1185 (ir->vdw_modifier != eintmodNONE && ir->rvdw == ir->rlist)) &&
1188 warning_note(wi, "With exact cut-offs, rlist should be "
1189 "larger than rcoulomb and rvdw, so that there "
1190 "is a buffer region for particle motion "
1191 "between neighborsearch steps");
1194 if (EEL_IS_ZERO_AT_CUTOFF(ir->coulombtype)
1195 && (ir->rlistlong <= ir->rcoulomb))
1197 sprintf(warn_buf, "For energy conservation with switch/shift potentials, %s should be 0.1 to 0.3 nm larger than rcoulomb.",
1198 IR_TWINRANGE(*ir) ? "rlistlong" : "rlist");
1199 warning_note(wi, warn_buf);
1201 if (EVDW_SWITCHED(ir->vdwtype) && (ir->rlistlong <= ir->rvdw))
1203 sprintf(warn_buf, "For energy conservation with switch/shift potentials, %s should be 0.1 to 0.3 nm larger than rvdw.",
1204 IR_TWINRANGE(*ir) ? "rlistlong" : "rlist");
1205 warning_note(wi, warn_buf);
1209 if (ir->vdwtype == evdwUSER && ir->eDispCorr != edispcNO)
1211 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.");
1214 if (ir->nstlist == -1)
1216 sprintf(err_buf, "With nstlist=-1 rvdw and rcoulomb should be smaller than rlist to account for diffusion and possibly charge-group radii");
1217 CHECK(ir->rvdw >= ir->rlist || ir->rcoulomb >= ir->rlist);
1219 sprintf(err_buf, "nstlist can not be smaller than -1");
1220 CHECK(ir->nstlist < -1);
1222 if (ir->eI == eiLBFGS && (ir->coulombtype == eelCUT || ir->vdwtype == evdwCUT)
1225 warning(wi, "For efficient BFGS minimization, use switch/shift/pme instead of cut-off.");
1228 if (ir->eI == eiLBFGS && ir->nbfgscorr <= 0)
1230 warning(wi, "Using L-BFGS with nbfgscorr<=0 just gets you steepest descent.");
1233 /* ENERGY CONSERVATION */
1234 if (ir_NVE(ir) && ir->cutoff_scheme == ecutsGROUP)
1236 if (!EVDW_MIGHT_BE_ZERO_AT_CUTOFF(ir->vdwtype) && ir->rvdw > 0 && ir->vdw_modifier == eintmodNONE)
1238 sprintf(warn_buf, "You are using a cut-off for VdW interactions with NVE, for good energy conservation use vdwtype = %s (possibly with DispCorr)",
1239 evdw_names[evdwSHIFT]);
1240 warning_note(wi, warn_buf);
1242 if (!EEL_MIGHT_BE_ZERO_AT_CUTOFF(ir->coulombtype) && ir->rcoulomb > 0 && ir->coulomb_modifier == eintmodNONE)
1244 sprintf(warn_buf, "You are using a cut-off for electrostatics with NVE, for good energy conservation use coulombtype = %s or %s",
1245 eel_names[eelPMESWITCH], eel_names[eelRF_ZERO]);
1246 warning_note(wi, warn_buf);
1250 /* IMPLICIT SOLVENT */
1251 if (ir->coulombtype == eelGB_NOTUSED)
1253 ir->coulombtype = eelCUT;
1254 ir->implicit_solvent = eisGBSA;
1255 fprintf(stderr, "Note: Old option for generalized born electrostatics given:\n"
1256 "Changing coulombtype from \"generalized-born\" to \"cut-off\" and instead\n"
1257 "setting implicit-solvent value to \"GBSA\" in input section.\n");
1260 if (ir->sa_algorithm == esaSTILL)
1262 sprintf(err_buf, "Still SA algorithm not available yet, use %s or %s instead\n", esa_names[esaAPPROX], esa_names[esaNO]);
1263 CHECK(ir->sa_algorithm == esaSTILL);
1266 if (ir->implicit_solvent == eisGBSA)
1268 sprintf(err_buf, "With GBSA implicit solvent, rgbradii must be equal to rlist.");
1269 CHECK(ir->rgbradii != ir->rlist);
1271 if (ir->coulombtype != eelCUT)
1273 sprintf(err_buf, "With GBSA, coulombtype must be equal to %s\n", eel_names[eelCUT]);
1274 CHECK(ir->coulombtype != eelCUT);
1276 if (ir->vdwtype != evdwCUT)
1278 sprintf(err_buf, "With GBSA, vdw-type must be equal to %s\n", evdw_names[evdwCUT]);
1279 CHECK(ir->vdwtype != evdwCUT);
1281 if (ir->nstgbradii < 1)
1283 sprintf(warn_buf, "Using GBSA with nstgbradii<1, setting nstgbradii=1");
1284 warning_note(wi, warn_buf);
1287 if (ir->sa_algorithm == esaNO)
1289 sprintf(warn_buf, "No SA (non-polar) calculation requested together with GB. Are you sure this is what you want?\n");
1290 warning_note(wi, warn_buf);
1292 if (ir->sa_surface_tension < 0 && ir->sa_algorithm != esaNO)
1294 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");
1295 warning_note(wi, warn_buf);
1297 if (ir->gb_algorithm == egbSTILL)
1299 ir->sa_surface_tension = 0.0049 * CAL2JOULE * 100;
1303 ir->sa_surface_tension = 0.0054 * CAL2JOULE * 100;
1306 if (ir->sa_surface_tension == 0 && ir->sa_algorithm != esaNO)
1308 sprintf(err_buf, "Surface tension set to 0 while SA-calculation requested\n");
1309 CHECK(ir->sa_surface_tension == 0 && ir->sa_algorithm != esaNO);
1316 if (ir->cutoff_scheme != ecutsGROUP)
1318 warning_error(wi, "AdresS simulation supports only cutoff-scheme=group");
1322 warning_error(wi, "AdresS simulation supports only stochastic dynamics");
1324 if (ir->epc != epcNO)
1326 warning_error(wi, "AdresS simulation does not support pressure coupling");
1328 if (EEL_FULL(ir->coulombtype))
1330 warning_error(wi, "AdresS simulation does not support long-range electrostatics");
1335 /* count the number of text elemets separated by whitespace in a string.
1336 str = the input string
1337 maxptr = the maximum number of allowed elements
1338 ptr = the output array of pointers to the first character of each element
1339 returns: the number of elements. */
1340 int str_nelem(const char *str, int maxptr, char *ptr[])
1345 copy0 = strdup(str);
1348 while (*copy != '\0')
1352 gmx_fatal(FARGS, "Too many groups on line: '%s' (max is %d)",
1360 while ((*copy != '\0') && !isspace(*copy))
1379 /* interpret a number of doubles from a string and put them in an array,
1380 after allocating space for them.
1381 str = the input string
1382 n = the (pre-allocated) number of doubles read
1383 r = the output array of doubles. */
1384 static void parse_n_real(char *str, int *n, real **r)
1389 *n = str_nelem(str, MAXPTR, ptr);
1392 for (i = 0; i < *n; i++)
1394 (*r)[i] = strtod(ptr[i], NULL);
1398 static void do_fep_params(t_inputrec *ir, char fep_lambda[][STRLEN], char weights[STRLEN])
1401 int i, j, max_n_lambda, nweights, nfep[efptNR];
1402 t_lambda *fep = ir->fepvals;
1403 t_expanded *expand = ir->expandedvals;
1404 real **count_fep_lambdas;
1405 gmx_bool bOneLambda = TRUE;
1407 snew(count_fep_lambdas, efptNR);
1409 /* FEP input processing */
1410 /* first, identify the number of lambda values for each type.
1411 All that are nonzero must have the same number */
1413 for (i = 0; i < efptNR; i++)
1415 parse_n_real(fep_lambda[i], &(nfep[i]), &(count_fep_lambdas[i]));
1418 /* now, determine the number of components. All must be either zero, or equal. */
1421 for (i = 0; i < efptNR; i++)
1423 if (nfep[i] > max_n_lambda)
1425 max_n_lambda = nfep[i]; /* here's a nonzero one. All of them
1426 must have the same number if its not zero.*/
1431 for (i = 0; i < efptNR; i++)
1435 ir->fepvals->separate_dvdl[i] = FALSE;
1437 else if (nfep[i] == max_n_lambda)
1439 if (i != efptTEMPERATURE) /* we treat this differently -- not really a reason to compute the derivative with
1440 respect to the temperature currently */
1442 ir->fepvals->separate_dvdl[i] = TRUE;
1447 gmx_fatal(FARGS, "Number of lambdas (%d) for FEP type %s not equal to number of other types (%d)",
1448 nfep[i], efpt_names[i], max_n_lambda);
1451 /* we don't print out dhdl if the temperature is changing, since we can't correctly define dhdl in this case */
1452 ir->fepvals->separate_dvdl[efptTEMPERATURE] = FALSE;
1454 /* the number of lambdas is the number we've read in, which is either zero
1455 or the same for all */
1456 fep->n_lambda = max_n_lambda;
1458 /* allocate space for the array of lambda values */
1459 snew(fep->all_lambda, efptNR);
1460 /* if init_lambda is defined, we need to set lambda */
1461 if ((fep->init_lambda > 0) && (fep->n_lambda == 0))
1463 ir->fepvals->separate_dvdl[efptFEP] = TRUE;
1465 /* otherwise allocate the space for all of the lambdas, and transfer the data */
1466 for (i = 0; i < efptNR; i++)
1468 snew(fep->all_lambda[i], fep->n_lambda);
1469 if (nfep[i] > 0) /* if it's zero, then the count_fep_lambda arrays
1472 for (j = 0; j < fep->n_lambda; j++)
1474 fep->all_lambda[i][j] = (double)count_fep_lambdas[i][j];
1476 sfree(count_fep_lambdas[i]);
1479 sfree(count_fep_lambdas);
1481 /* "fep-vals" is either zero or the full number. If zero, we'll need to define fep-lambdas for internal
1482 bookkeeping -- for now, init_lambda */
1484 if ((nfep[efptFEP] == 0) && (fep->init_lambda >= 0))
1486 for (i = 0; i < fep->n_lambda; i++)
1488 fep->all_lambda[efptFEP][i] = fep->init_lambda;
1492 /* check to see if only a single component lambda is defined, and soft core is defined.
1493 In this case, turn on coulomb soft core */
1495 if (max_n_lambda == 0)
1501 for (i = 0; i < efptNR; i++)
1503 if ((nfep[i] != 0) && (i != efptFEP))
1509 if ((bOneLambda) && (fep->sc_alpha > 0))
1511 fep->bScCoul = TRUE;
1514 /* Fill in the others with the efptFEP if they are not explicitly
1515 specified (i.e. nfep[i] == 0). This means if fep is not defined,
1516 they are all zero. */
1518 for (i = 0; i < efptNR; i++)
1520 if ((nfep[i] == 0) && (i != efptFEP))
1522 for (j = 0; j < fep->n_lambda; j++)
1524 fep->all_lambda[i][j] = fep->all_lambda[efptFEP][j];
1530 /* make it easier if sc_r_power = 48 by increasing it to the 4th power, to be in the right scale. */
1531 if (fep->sc_r_power == 48)
1533 if (fep->sc_alpha > 0.1)
1535 gmx_fatal(FARGS, "sc_alpha (%f) for sc_r_power = 48 should usually be between 0.001 and 0.004", fep->sc_alpha);
1539 expand = ir->expandedvals;
1540 /* now read in the weights */
1541 parse_n_real(weights, &nweights, &(expand->init_lambda_weights));
1544 snew(expand->init_lambda_weights, fep->n_lambda); /* initialize to zero */
1546 else if (nweights != fep->n_lambda)
1548 gmx_fatal(FARGS, "Number of weights (%d) is not equal to number of lambda values (%d)",
1549 nweights, fep->n_lambda);
1551 if ((expand->nstexpanded < 0) && (ir->efep != efepNO))
1553 expand->nstexpanded = fep->nstdhdl;
1554 /* if you don't specify nstexpanded when doing expanded ensemble free energy calcs, it is set to nstdhdl */
1556 if ((expand->nstexpanded < 0) && ir->bSimTemp)
1558 expand->nstexpanded = 2*(int)(ir->opts.tau_t[0]/ir->delta_t);
1559 /* if you don't specify nstexpanded when doing expanded ensemble simulated tempering, it is set to
1560 2*tau_t just to be careful so it's not to frequent */
1565 static void do_simtemp_params(t_inputrec *ir)
1568 snew(ir->simtempvals->temperatures, ir->fepvals->n_lambda);
1569 GetSimTemps(ir->fepvals->n_lambda, ir->simtempvals, ir->fepvals->all_lambda[efptTEMPERATURE]);
1574 static void do_wall_params(t_inputrec *ir,
1575 char *wall_atomtype, char *wall_density,
1579 char *names[MAXPTR];
1582 opts->wall_atomtype[0] = NULL;
1583 opts->wall_atomtype[1] = NULL;
1585 ir->wall_atomtype[0] = -1;
1586 ir->wall_atomtype[1] = -1;
1587 ir->wall_density[0] = 0;
1588 ir->wall_density[1] = 0;
1592 nstr = str_nelem(wall_atomtype, MAXPTR, names);
1593 if (nstr != ir->nwall)
1595 gmx_fatal(FARGS, "Expected %d elements for wall_atomtype, found %d",
1598 for (i = 0; i < ir->nwall; i++)
1600 opts->wall_atomtype[i] = strdup(names[i]);
1603 if (ir->wall_type == ewt93 || ir->wall_type == ewt104)
1605 nstr = str_nelem(wall_density, MAXPTR, names);
1606 if (nstr != ir->nwall)
1608 gmx_fatal(FARGS, "Expected %d elements for wall-density, found %d", ir->nwall, nstr);
1610 for (i = 0; i < ir->nwall; i++)
1612 sscanf(names[i], "%lf", &dbl);
1615 gmx_fatal(FARGS, "wall-density[%d] = %f\n", i, dbl);
1617 ir->wall_density[i] = dbl;
1623 static void add_wall_energrps(gmx_groups_t *groups, int nwall, t_symtab *symtab)
1631 srenew(groups->grpname, groups->ngrpname+nwall);
1632 grps = &(groups->grps[egcENER]);
1633 srenew(grps->nm_ind, grps->nr+nwall);
1634 for (i = 0; i < nwall; i++)
1636 sprintf(str, "wall%d", i);
1637 groups->grpname[groups->ngrpname] = put_symtab(symtab, str);
1638 grps->nm_ind[grps->nr++] = groups->ngrpname++;
1643 void read_expandedparams(int *ninp_p, t_inpfile **inp_p,
1644 t_expanded *expand, warninp_t wi)
1646 int ninp, nerror = 0;
1652 /* read expanded ensemble parameters */
1653 CCTYPE ("expanded ensemble variables");
1654 ITYPE ("nstexpanded", expand->nstexpanded, -1);
1655 EETYPE("lmc-stats", expand->elamstats, elamstats_names);
1656 EETYPE("lmc-move", expand->elmcmove, elmcmove_names);
1657 EETYPE("lmc-weights-equil", expand->elmceq, elmceq_names);
1658 ITYPE ("weight-equil-number-all-lambda", expand->equil_n_at_lam, -1);
1659 ITYPE ("weight-equil-number-samples", expand->equil_samples, -1);
1660 ITYPE ("weight-equil-number-steps", expand->equil_steps, -1);
1661 RTYPE ("weight-equil-wl-delta", expand->equil_wl_delta, -1);
1662 RTYPE ("weight-equil-count-ratio", expand->equil_ratio, -1);
1663 CCTYPE("Seed for Monte Carlo in lambda space");
1664 ITYPE ("lmc-seed", expand->lmc_seed, -1);
1665 RTYPE ("mc-temperature", expand->mc_temp, -1);
1666 ITYPE ("lmc-repeats", expand->lmc_repeats, 1);
1667 ITYPE ("lmc-gibbsdelta", expand->gibbsdeltalam, -1);
1668 ITYPE ("lmc-forced-nstart", expand->lmc_forced_nstart, 0);
1669 EETYPE("symmetrized-transition-matrix", expand->bSymmetrizedTMatrix, yesno_names);
1670 ITYPE("nst-transition-matrix", expand->nstTij, -1);
1671 ITYPE ("mininum-var-min", expand->minvarmin, 100); /*default is reasonable */
1672 ITYPE ("weight-c-range", expand->c_range, 0); /* default is just C=0 */
1673 RTYPE ("wl-scale", expand->wl_scale, 0.8);
1674 RTYPE ("wl-ratio", expand->wl_ratio, 0.8);
1675 RTYPE ("init-wl-delta", expand->init_wl_delta, 1.0);
1676 EETYPE("wl-oneovert", expand->bWLoneovert, yesno_names);
1684 void get_ir(const char *mdparin, const char *mdparout,
1685 t_inputrec *ir, t_gromppopts *opts,
1689 double dumdub[2][6];
1693 char warn_buf[STRLEN];
1694 t_lambda *fep = ir->fepvals;
1695 t_expanded *expand = ir->expandedvals;
1697 inp = read_inpfile(mdparin, &ninp, wi);
1699 snew(dumstr[0], STRLEN);
1700 snew(dumstr[1], STRLEN);
1702 if (-1 == search_einp(ninp, inp, "cutoff-scheme"))
1705 "%s did not specify a value for the .mdp option "
1706 "\"cutoff-scheme\". Probably it was first intended for use "
1707 "with GROMACS before 4.6. In 4.6, the Verlet scheme was "
1708 "introduced, but the group scheme was still the default. "
1709 "The default is now the Verlet scheme, so you will observe "
1710 "different behaviour.", mdparin);
1711 warning_note(wi, warn_buf);
1714 /* remove the following deprecated commands */
1717 REM_TYPE("domain-decomposition");
1718 REM_TYPE("andersen-seed");
1720 REM_TYPE("dihre-fc");
1721 REM_TYPE("dihre-tau");
1722 REM_TYPE("nstdihreout");
1723 REM_TYPE("nstcheckpoint");
1725 /* replace the following commands with the clearer new versions*/
1726 REPL_TYPE("unconstrained-start", "continuation");
1727 REPL_TYPE("foreign-lambda", "fep-lambdas");
1728 REPL_TYPE("verlet-buffer-drift", "verlet-buffer-tolerance");
1730 CCTYPE ("VARIOUS PREPROCESSING OPTIONS");
1731 CTYPE ("Preprocessor information: use cpp syntax.");
1732 CTYPE ("e.g.: -I/home/joe/doe -I/home/mary/roe");
1733 STYPE ("include", opts->include, NULL);
1734 CTYPE ("e.g.: -DPOSRES -DFLEXIBLE (note these variable names are case sensitive)");
1735 STYPE ("define", opts->define, NULL);
1737 CCTYPE ("RUN CONTROL PARAMETERS");
1738 EETYPE("integrator", ir->eI, ei_names);
1739 CTYPE ("Start time and timestep in ps");
1740 RTYPE ("tinit", ir->init_t, 0.0);
1741 RTYPE ("dt", ir->delta_t, 0.001);
1742 STEPTYPE ("nsteps", ir->nsteps, 0);
1743 CTYPE ("For exact run continuation or redoing part of a run");
1744 STEPTYPE ("init-step", ir->init_step, 0);
1745 CTYPE ("Part index is updated automatically on checkpointing (keeps files separate)");
1746 ITYPE ("simulation-part", ir->simulation_part, 1);
1747 CTYPE ("mode for center of mass motion removal");
1748 EETYPE("comm-mode", ir->comm_mode, ecm_names);
1749 CTYPE ("number of steps for center of mass motion removal");
1750 ITYPE ("nstcomm", ir->nstcomm, 100);
1751 CTYPE ("group(s) for center of mass motion removal");
1752 STYPE ("comm-grps", vcm, NULL);
1754 CCTYPE ("LANGEVIN DYNAMICS OPTIONS");
1755 CTYPE ("Friction coefficient (amu/ps) and random seed");
1756 RTYPE ("bd-fric", ir->bd_fric, 0.0);
1757 ITYPE ("ld-seed", ir->ld_seed, 1993);
1760 CCTYPE ("ENERGY MINIMIZATION OPTIONS");
1761 CTYPE ("Force tolerance and initial step-size");
1762 RTYPE ("emtol", ir->em_tol, 10.0);
1763 RTYPE ("emstep", ir->em_stepsize, 0.01);
1764 CTYPE ("Max number of iterations in relax-shells");
1765 ITYPE ("niter", ir->niter, 20);
1766 CTYPE ("Step size (ps^2) for minimization of flexible constraints");
1767 RTYPE ("fcstep", ir->fc_stepsize, 0);
1768 CTYPE ("Frequency of steepest descents steps when doing CG");
1769 ITYPE ("nstcgsteep", ir->nstcgsteep, 1000);
1770 ITYPE ("nbfgscorr", ir->nbfgscorr, 10);
1772 CCTYPE ("TEST PARTICLE INSERTION OPTIONS");
1773 RTYPE ("rtpi", ir->rtpi, 0.05);
1775 /* Output options */
1776 CCTYPE ("OUTPUT CONTROL OPTIONS");
1777 CTYPE ("Output frequency for coords (x), velocities (v) and forces (f)");
1778 ITYPE ("nstxout", ir->nstxout, 0);
1779 ITYPE ("nstvout", ir->nstvout, 0);
1780 ITYPE ("nstfout", ir->nstfout, 0);
1781 ir->nstcheckpoint = 1000;
1782 CTYPE ("Output frequency for energies to log file and energy file");
1783 ITYPE ("nstlog", ir->nstlog, 1000);
1784 ITYPE ("nstcalcenergy", ir->nstcalcenergy, 100);
1785 ITYPE ("nstenergy", ir->nstenergy, 1000);
1786 CTYPE ("Output frequency and precision for .xtc file");
1787 ITYPE ("nstxtcout", ir->nstxtcout, 0);
1788 RTYPE ("xtc-precision", ir->xtcprec, 1000.0);
1789 CTYPE ("This selects the subset of atoms for the .xtc file. You can");
1790 CTYPE ("select multiple groups. By default all atoms will be written.");
1791 STYPE ("xtc-grps", xtc_grps, NULL);
1792 CTYPE ("Selection of energy groups");
1793 STYPE ("energygrps", energy, NULL);
1795 /* Neighbor searching */
1796 CCTYPE ("NEIGHBORSEARCHING PARAMETERS");
1797 CTYPE ("cut-off scheme (Verlet: particle based cut-offs, group: using charge groups)");
1798 EETYPE("cutoff-scheme", ir->cutoff_scheme, ecutscheme_names);
1799 CTYPE ("nblist update frequency");
1800 ITYPE ("nstlist", ir->nstlist, 10);
1801 CTYPE ("ns algorithm (simple or grid)");
1802 EETYPE("ns-type", ir->ns_type, ens_names);
1803 /* set ndelta to the optimal value of 2 */
1805 CTYPE ("Periodic boundary conditions: xyz, no, xy");
1806 EETYPE("pbc", ir->ePBC, epbc_names);
1807 EETYPE("periodic-molecules", ir->bPeriodicMols, yesno_names);
1808 CTYPE ("Allowed energy error due to the Verlet buffer in kJ/mol/ps per atom,");
1809 CTYPE ("a value of -1 means: use rlist");
1810 RTYPE("verlet-buffer-tolerance", ir->verletbuf_tol, 0.005);
1811 CTYPE ("nblist cut-off");
1812 RTYPE ("rlist", ir->rlist, 1.0);
1813 CTYPE ("long-range cut-off for switched potentials");
1814 RTYPE ("rlistlong", ir->rlistlong, -1);
1815 ITYPE ("nstcalclr", ir->nstcalclr, -1);
1817 /* Electrostatics */
1818 CCTYPE ("OPTIONS FOR ELECTROSTATICS AND VDW");
1819 CTYPE ("Method for doing electrostatics");
1820 EETYPE("coulombtype", ir->coulombtype, eel_names);
1821 EETYPE("coulomb-modifier", ir->coulomb_modifier, eintmod_names);
1822 CTYPE ("cut-off lengths");
1823 RTYPE ("rcoulomb-switch", ir->rcoulomb_switch, 0.0);
1824 RTYPE ("rcoulomb", ir->rcoulomb, 1.0);
1825 CTYPE ("Relative dielectric constant for the medium and the reaction field");
1826 RTYPE ("epsilon-r", ir->epsilon_r, 1.0);
1827 RTYPE ("epsilon-rf", ir->epsilon_rf, 0.0);
1828 CTYPE ("Method for doing Van der Waals");
1829 EETYPE("vdw-type", ir->vdwtype, evdw_names);
1830 EETYPE("vdw-modifier", ir->vdw_modifier, eintmod_names);
1831 CTYPE ("cut-off lengths");
1832 RTYPE ("rvdw-switch", ir->rvdw_switch, 0.0);
1833 RTYPE ("rvdw", ir->rvdw, 1.0);
1834 CTYPE ("Apply long range dispersion corrections for Energy and Pressure");
1835 EETYPE("DispCorr", ir->eDispCorr, edispc_names);
1836 CTYPE ("Extension of the potential lookup tables beyond the cut-off");
1837 RTYPE ("table-extension", ir->tabext, 1.0);
1838 CTYPE ("Separate tables between energy group pairs");
1839 STYPE ("energygrp-table", egptable, NULL);
1840 CTYPE ("Spacing for the PME/PPPM FFT grid");
1841 RTYPE ("fourierspacing", ir->fourier_spacing, 0.12);
1842 CTYPE ("FFT grid size, when a value is 0 fourierspacing will be used");
1843 ITYPE ("fourier-nx", ir->nkx, 0);
1844 ITYPE ("fourier-ny", ir->nky, 0);
1845 ITYPE ("fourier-nz", ir->nkz, 0);
1846 CTYPE ("EWALD/PME/PPPM parameters");
1847 ITYPE ("pme-order", ir->pme_order, 4);
1848 RTYPE ("ewald-rtol", ir->ewald_rtol, 0.00001);
1849 RTYPE ("ewald-rtol-lj", ir->ewald_rtol_lj, 0.001);
1850 EETYPE("lj-pme-comb-rule", ir->ljpme_combination_rule, eljpme_names);
1851 EETYPE("ewald-geometry", ir->ewald_geometry, eewg_names);
1852 RTYPE ("epsilon-surface", ir->epsilon_surface, 0.0);
1853 EETYPE("optimize-fft", ir->bOptFFT, yesno_names);
1855 CCTYPE("IMPLICIT SOLVENT ALGORITHM");
1856 EETYPE("implicit-solvent", ir->implicit_solvent, eis_names);
1858 CCTYPE ("GENERALIZED BORN ELECTROSTATICS");
1859 CTYPE ("Algorithm for calculating Born radii");
1860 EETYPE("gb-algorithm", ir->gb_algorithm, egb_names);
1861 CTYPE ("Frequency of calculating the Born radii inside rlist");
1862 ITYPE ("nstgbradii", ir->nstgbradii, 1);
1863 CTYPE ("Cutoff for Born radii calculation; the contribution from atoms");
1864 CTYPE ("between rlist and rgbradii is updated every nstlist steps");
1865 RTYPE ("rgbradii", ir->rgbradii, 1.0);
1866 CTYPE ("Dielectric coefficient of the implicit solvent");
1867 RTYPE ("gb-epsilon-solvent", ir->gb_epsilon_solvent, 80.0);
1868 CTYPE ("Salt concentration in M for Generalized Born models");
1869 RTYPE ("gb-saltconc", ir->gb_saltconc, 0.0);
1870 CTYPE ("Scaling factors used in the OBC GB model. Default values are OBC(II)");
1871 RTYPE ("gb-obc-alpha", ir->gb_obc_alpha, 1.0);
1872 RTYPE ("gb-obc-beta", ir->gb_obc_beta, 0.8);
1873 RTYPE ("gb-obc-gamma", ir->gb_obc_gamma, 4.85);
1874 RTYPE ("gb-dielectric-offset", ir->gb_dielectric_offset, 0.009);
1875 EETYPE("sa-algorithm", ir->sa_algorithm, esa_names);
1876 CTYPE ("Surface tension (kJ/mol/nm^2) for the SA (nonpolar surface) part of GBSA");
1877 CTYPE ("The value -1 will set default value for Still/HCT/OBC GB-models.");
1878 RTYPE ("sa-surface-tension", ir->sa_surface_tension, -1);
1880 /* Coupling stuff */
1881 CCTYPE ("OPTIONS FOR WEAK COUPLING ALGORITHMS");
1882 CTYPE ("Temperature coupling");
1883 EETYPE("tcoupl", ir->etc, etcoupl_names);
1884 ITYPE ("nsttcouple", ir->nsttcouple, -1);
1885 ITYPE("nh-chain-length", ir->opts.nhchainlength, NHCHAINLENGTH);
1886 EETYPE("print-nose-hoover-chain-variables", ir->bPrintNHChains, yesno_names);
1887 CTYPE ("Groups to couple separately");
1888 STYPE ("tc-grps", tcgrps, NULL);
1889 CTYPE ("Time constant (ps) and reference temperature (K)");
1890 STYPE ("tau-t", tau_t, NULL);
1891 STYPE ("ref-t", ref_t, NULL);
1892 CTYPE ("pressure coupling");
1893 EETYPE("pcoupl", ir->epc, epcoupl_names);
1894 EETYPE("pcoupltype", ir->epct, epcoupltype_names);
1895 ITYPE ("nstpcouple", ir->nstpcouple, -1);
1896 CTYPE ("Time constant (ps), compressibility (1/bar) and reference P (bar)");
1897 RTYPE ("tau-p", ir->tau_p, 1.0);
1898 STYPE ("compressibility", dumstr[0], NULL);
1899 STYPE ("ref-p", dumstr[1], NULL);
1900 CTYPE ("Scaling of reference coordinates, No, All or COM");
1901 EETYPE ("refcoord-scaling", ir->refcoord_scaling, erefscaling_names);
1904 CCTYPE ("OPTIONS FOR QMMM calculations");
1905 EETYPE("QMMM", ir->bQMMM, yesno_names);
1906 CTYPE ("Groups treated Quantum Mechanically");
1907 STYPE ("QMMM-grps", QMMM, NULL);
1908 CTYPE ("QM method");
1909 STYPE("QMmethod", QMmethod, NULL);
1910 CTYPE ("QMMM scheme");
1911 EETYPE("QMMMscheme", ir->QMMMscheme, eQMMMscheme_names);
1912 CTYPE ("QM basisset");
1913 STYPE("QMbasis", QMbasis, NULL);
1914 CTYPE ("QM charge");
1915 STYPE ("QMcharge", QMcharge, NULL);
1916 CTYPE ("QM multiplicity");
1917 STYPE ("QMmult", QMmult, NULL);
1918 CTYPE ("Surface Hopping");
1919 STYPE ("SH", bSH, NULL);
1920 CTYPE ("CAS space options");
1921 STYPE ("CASorbitals", CASorbitals, NULL);
1922 STYPE ("CASelectrons", CASelectrons, NULL);
1923 STYPE ("SAon", SAon, NULL);
1924 STYPE ("SAoff", SAoff, NULL);
1925 STYPE ("SAsteps", SAsteps, NULL);
1926 CTYPE ("Scale factor for MM charges");
1927 RTYPE ("MMChargeScaleFactor", ir->scalefactor, 1.0);
1928 CTYPE ("Optimization of QM subsystem");
1929 STYPE ("bOPT", bOPT, NULL);
1930 STYPE ("bTS", bTS, NULL);
1932 /* Simulated annealing */
1933 CCTYPE("SIMULATED ANNEALING");
1934 CTYPE ("Type of annealing for each temperature group (no/single/periodic)");
1935 STYPE ("annealing", anneal, NULL);
1936 CTYPE ("Number of time points to use for specifying annealing in each group");
1937 STYPE ("annealing-npoints", anneal_npoints, NULL);
1938 CTYPE ("List of times at the annealing points for each group");
1939 STYPE ("annealing-time", anneal_time, NULL);
1940 CTYPE ("Temp. at each annealing point, for each group.");
1941 STYPE ("annealing-temp", anneal_temp, NULL);
1944 CCTYPE ("GENERATE VELOCITIES FOR STARTUP RUN");
1945 EETYPE("gen-vel", opts->bGenVel, yesno_names);
1946 RTYPE ("gen-temp", opts->tempi, 300.0);
1947 ITYPE ("gen-seed", opts->seed, 173529);
1950 CCTYPE ("OPTIONS FOR BONDS");
1951 EETYPE("constraints", opts->nshake, constraints);
1952 CTYPE ("Type of constraint algorithm");
1953 EETYPE("constraint-algorithm", ir->eConstrAlg, econstr_names);
1954 CTYPE ("Do not constrain the start configuration");
1955 EETYPE("continuation", ir->bContinuation, yesno_names);
1956 CTYPE ("Use successive overrelaxation to reduce the number of shake iterations");
1957 EETYPE("Shake-SOR", ir->bShakeSOR, yesno_names);
1958 CTYPE ("Relative tolerance of shake");
1959 RTYPE ("shake-tol", ir->shake_tol, 0.0001);
1960 CTYPE ("Highest order in the expansion of the constraint coupling matrix");
1961 ITYPE ("lincs-order", ir->nProjOrder, 4);
1962 CTYPE ("Number of iterations in the final step of LINCS. 1 is fine for");
1963 CTYPE ("normal simulations, but use 2 to conserve energy in NVE runs.");
1964 CTYPE ("For energy minimization with constraints it should be 4 to 8.");
1965 ITYPE ("lincs-iter", ir->nLincsIter, 1);
1966 CTYPE ("Lincs will write a warning to the stderr if in one step a bond");
1967 CTYPE ("rotates over more degrees than");
1968 RTYPE ("lincs-warnangle", ir->LincsWarnAngle, 30.0);
1969 CTYPE ("Convert harmonic bonds to morse potentials");
1970 EETYPE("morse", opts->bMorse, yesno_names);
1972 /* Energy group exclusions */
1973 CCTYPE ("ENERGY GROUP EXCLUSIONS");
1974 CTYPE ("Pairs of energy groups for which all non-bonded interactions are excluded");
1975 STYPE ("energygrp-excl", egpexcl, NULL);
1979 CTYPE ("Number of walls, type, atom types, densities and box-z scale factor for Ewald");
1980 ITYPE ("nwall", ir->nwall, 0);
1981 EETYPE("wall-type", ir->wall_type, ewt_names);
1982 RTYPE ("wall-r-linpot", ir->wall_r_linpot, -1);
1983 STYPE ("wall-atomtype", wall_atomtype, NULL);
1984 STYPE ("wall-density", wall_density, NULL);
1985 RTYPE ("wall-ewald-zfac", ir->wall_ewald_zfac, 3);
1988 CCTYPE("COM PULLING");
1989 CTYPE("Pull type: no, umbrella, constraint or constant-force");
1990 EETYPE("pull", ir->ePull, epull_names);
1991 if (ir->ePull != epullNO)
1994 pull_grp = read_pullparams(&ninp, &inp, ir->pull, &opts->pull_start, wi);
1997 /* Enforced rotation */
1998 CCTYPE("ENFORCED ROTATION");
1999 CTYPE("Enforced rotation: No or Yes");
2000 EETYPE("rotation", ir->bRot, yesno_names);
2004 rot_grp = read_rotparams(&ninp, &inp, ir->rot, wi);
2008 CCTYPE("NMR refinement stuff");
2009 CTYPE ("Distance restraints type: No, Simple or Ensemble");
2010 EETYPE("disre", ir->eDisre, edisre_names);
2011 CTYPE ("Force weighting of pairs in one distance restraint: Conservative or Equal");
2012 EETYPE("disre-weighting", ir->eDisreWeighting, edisreweighting_names);
2013 CTYPE ("Use sqrt of the time averaged times the instantaneous violation");
2014 EETYPE("disre-mixed", ir->bDisreMixed, yesno_names);
2015 RTYPE ("disre-fc", ir->dr_fc, 1000.0);
2016 RTYPE ("disre-tau", ir->dr_tau, 0.0);
2017 CTYPE ("Output frequency for pair distances to energy file");
2018 ITYPE ("nstdisreout", ir->nstdisreout, 100);
2019 CTYPE ("Orientation restraints: No or Yes");
2020 EETYPE("orire", opts->bOrire, yesno_names);
2021 CTYPE ("Orientation restraints force constant and tau for time averaging");
2022 RTYPE ("orire-fc", ir->orires_fc, 0.0);
2023 RTYPE ("orire-tau", ir->orires_tau, 0.0);
2024 STYPE ("orire-fitgrp", orirefitgrp, NULL);
2025 CTYPE ("Output frequency for trace(SD) and S to energy file");
2026 ITYPE ("nstorireout", ir->nstorireout, 100);
2028 /* free energy variables */
2029 CCTYPE ("Free energy variables");
2030 EETYPE("free-energy", ir->efep, efep_names);
2031 STYPE ("couple-moltype", couple_moltype, NULL);
2032 EETYPE("couple-lambda0", opts->couple_lam0, couple_lam);
2033 EETYPE("couple-lambda1", opts->couple_lam1, couple_lam);
2034 EETYPE("couple-intramol", opts->bCoupleIntra, yesno_names);
2036 RTYPE ("init-lambda", fep->init_lambda, -1); /* start with -1 so
2038 it was not entered */
2039 ITYPE ("init-lambda-state", fep->init_fep_state, -1);
2040 RTYPE ("delta-lambda", fep->delta_lambda, 0.0);
2041 ITYPE ("nstdhdl", fep->nstdhdl, 50);
2042 STYPE ("fep-lambdas", fep_lambda[efptFEP], NULL);
2043 STYPE ("mass-lambdas", fep_lambda[efptMASS], NULL);
2044 STYPE ("coul-lambdas", fep_lambda[efptCOUL], NULL);
2045 STYPE ("vdw-lambdas", fep_lambda[efptVDW], NULL);
2046 STYPE ("bonded-lambdas", fep_lambda[efptBONDED], NULL);
2047 STYPE ("restraint-lambdas", fep_lambda[efptRESTRAINT], NULL);
2048 STYPE ("temperature-lambdas", fep_lambda[efptTEMPERATURE], NULL);
2049 ITYPE ("calc-lambda-neighbors", fep->lambda_neighbors, 1);
2050 STYPE ("init-lambda-weights", lambda_weights, NULL);
2051 EETYPE("dhdl-print-energy", fep->bPrintEnergy, yesno_names);
2052 RTYPE ("sc-alpha", fep->sc_alpha, 0.0);
2053 ITYPE ("sc-power", fep->sc_power, 1);
2054 RTYPE ("sc-r-power", fep->sc_r_power, 6.0);
2055 RTYPE ("sc-sigma", fep->sc_sigma, 0.3);
2056 EETYPE("sc-coul", fep->bScCoul, yesno_names);
2057 ITYPE ("dh_hist_size", fep->dh_hist_size, 0);
2058 RTYPE ("dh_hist_spacing", fep->dh_hist_spacing, 0.1);
2059 EETYPE("separate-dhdl-file", fep->separate_dhdl_file,
2060 separate_dhdl_file_names);
2061 EETYPE("dhdl-derivatives", fep->dhdl_derivatives, dhdl_derivatives_names);
2062 ITYPE ("dh_hist_size", fep->dh_hist_size, 0);
2063 RTYPE ("dh_hist_spacing", fep->dh_hist_spacing, 0.1);
2065 /* Non-equilibrium MD stuff */
2066 CCTYPE("Non-equilibrium MD stuff");
2067 STYPE ("acc-grps", accgrps, NULL);
2068 STYPE ("accelerate", acc, NULL);
2069 STYPE ("freezegrps", freeze, NULL);
2070 STYPE ("freezedim", frdim, NULL);
2071 RTYPE ("cos-acceleration", ir->cos_accel, 0);
2072 STYPE ("deform", deform, NULL);
2074 /* simulated tempering variables */
2075 CCTYPE("simulated tempering variables");
2076 EETYPE("simulated-tempering", ir->bSimTemp, yesno_names);
2077 EETYPE("simulated-tempering-scaling", ir->simtempvals->eSimTempScale, esimtemp_names);
2078 RTYPE("sim-temp-low", ir->simtempvals->simtemp_low, 300.0);
2079 RTYPE("sim-temp-high", ir->simtempvals->simtemp_high, 300.0);
2081 /* expanded ensemble variables */
2082 if (ir->efep == efepEXPANDED || ir->bSimTemp)
2084 read_expandedparams(&ninp, &inp, expand, wi);
2087 /* Electric fields */
2088 CCTYPE("Electric fields");
2089 CTYPE ("Format is number of terms (int) and for all terms an amplitude (real)");
2090 CTYPE ("and a phase angle (real)");
2091 STYPE ("E-x", efield_x, NULL);
2092 STYPE ("E-xt", efield_xt, NULL);
2093 STYPE ("E-y", efield_y, NULL);
2094 STYPE ("E-yt", efield_yt, NULL);
2095 STYPE ("E-z", efield_z, NULL);
2096 STYPE ("E-zt", efield_zt, NULL);
2098 /* AdResS defined thingies */
2099 CCTYPE ("AdResS parameters");
2100 EETYPE("adress", ir->bAdress, yesno_names);
2103 snew(ir->adress, 1);
2104 read_adressparams(&ninp, &inp, ir->adress, wi);
2107 /* User defined thingies */
2108 CCTYPE ("User defined thingies");
2109 STYPE ("user1-grps", user1, NULL);
2110 STYPE ("user2-grps", user2, NULL);
2111 ITYPE ("userint1", ir->userint1, 0);
2112 ITYPE ("userint2", ir->userint2, 0);
2113 ITYPE ("userint3", ir->userint3, 0);
2114 ITYPE ("userint4", ir->userint4, 0);
2115 RTYPE ("userreal1", ir->userreal1, 0);
2116 RTYPE ("userreal2", ir->userreal2, 0);
2117 RTYPE ("userreal3", ir->userreal3, 0);
2118 RTYPE ("userreal4", ir->userreal4, 0);
2121 write_inpfile(mdparout, ninp, inp, FALSE, wi);
2122 for (i = 0; (i < ninp); i++)
2125 sfree(inp[i].value);
2129 /* Process options if necessary */
2130 for (m = 0; m < 2; m++)
2132 for (i = 0; i < 2*DIM; i++)
2141 if (sscanf(dumstr[m], "%lf", &(dumdub[m][XX])) != 1)
2143 warning_error(wi, "Pressure coupling not enough values (I need 1)");
2145 dumdub[m][YY] = dumdub[m][ZZ] = dumdub[m][XX];
2147 case epctSEMIISOTROPIC:
2148 case epctSURFACETENSION:
2149 if (sscanf(dumstr[m], "%lf%lf",
2150 &(dumdub[m][XX]), &(dumdub[m][ZZ])) != 2)
2152 warning_error(wi, "Pressure coupling not enough values (I need 2)");
2154 dumdub[m][YY] = dumdub[m][XX];
2156 case epctANISOTROPIC:
2157 if (sscanf(dumstr[m], "%lf%lf%lf%lf%lf%lf",
2158 &(dumdub[m][XX]), &(dumdub[m][YY]), &(dumdub[m][ZZ]),
2159 &(dumdub[m][3]), &(dumdub[m][4]), &(dumdub[m][5])) != 6)
2161 warning_error(wi, "Pressure coupling not enough values (I need 6)");
2165 gmx_fatal(FARGS, "Pressure coupling type %s not implemented yet",
2166 epcoupltype_names[ir->epct]);
2170 clear_mat(ir->ref_p);
2171 clear_mat(ir->compress);
2172 for (i = 0; i < DIM; i++)
2174 ir->ref_p[i][i] = dumdub[1][i];
2175 ir->compress[i][i] = dumdub[0][i];
2177 if (ir->epct == epctANISOTROPIC)
2179 ir->ref_p[XX][YY] = dumdub[1][3];
2180 ir->ref_p[XX][ZZ] = dumdub[1][4];
2181 ir->ref_p[YY][ZZ] = dumdub[1][5];
2182 if (ir->ref_p[XX][YY] != 0 && ir->ref_p[XX][ZZ] != 0 && ir->ref_p[YY][ZZ] != 0)
2184 warning(wi, "All off-diagonal reference pressures are non-zero. Are you sure you want to apply a threefold shear stress?\n");
2186 ir->compress[XX][YY] = dumdub[0][3];
2187 ir->compress[XX][ZZ] = dumdub[0][4];
2188 ir->compress[YY][ZZ] = dumdub[0][5];
2189 for (i = 0; i < DIM; i++)
2191 for (m = 0; m < i; m++)
2193 ir->ref_p[i][m] = ir->ref_p[m][i];
2194 ir->compress[i][m] = ir->compress[m][i];
2199 if (ir->comm_mode == ecmNO)
2204 opts->couple_moltype = NULL;
2205 if (strlen(couple_moltype) > 0)
2207 if (ir->efep != efepNO)
2209 opts->couple_moltype = strdup(couple_moltype);
2210 if (opts->couple_lam0 == opts->couple_lam1)
2212 warning(wi, "The lambda=0 and lambda=1 states for coupling are identical");
2214 if (ir->eI == eiMD && (opts->couple_lam0 == ecouplamNONE ||
2215 opts->couple_lam1 == ecouplamNONE))
2217 warning(wi, "For proper sampling of the (nearly) decoupled state, stochastic dynamics should be used");
2222 warning(wi, "Can not couple a molecule with free_energy = no");
2225 /* FREE ENERGY AND EXPANDED ENSEMBLE OPTIONS */
2226 if (ir->efep != efepNO)
2228 if (fep->delta_lambda > 0)
2230 ir->efep = efepSLOWGROWTH;
2236 fep->bPrintEnergy = TRUE;
2237 /* always print out the energy to dhdl if we are doing expanded ensemble, since we need the total energy
2238 if the temperature is changing. */
2241 if ((ir->efep != efepNO) || ir->bSimTemp)
2243 ir->bExpanded = FALSE;
2244 if ((ir->efep == efepEXPANDED) || ir->bSimTemp)
2246 ir->bExpanded = TRUE;
2248 do_fep_params(ir, fep_lambda, lambda_weights);
2249 if (ir->bSimTemp) /* done after fep params */
2251 do_simtemp_params(ir);
2256 ir->fepvals->n_lambda = 0;
2259 /* WALL PARAMETERS */
2261 do_wall_params(ir, wall_atomtype, wall_density, opts);
2263 /* ORIENTATION RESTRAINT PARAMETERS */
2265 if (opts->bOrire && str_nelem(orirefitgrp, MAXPTR, NULL) != 1)
2267 warning_error(wi, "ERROR: Need one orientation restraint fit group\n");
2270 /* DEFORMATION PARAMETERS */
2272 clear_mat(ir->deform);
2273 for (i = 0; i < 6; i++)
2277 m = sscanf(deform, "%lf %lf %lf %lf %lf %lf",
2278 &(dumdub[0][0]), &(dumdub[0][1]), &(dumdub[0][2]),
2279 &(dumdub[0][3]), &(dumdub[0][4]), &(dumdub[0][5]));
2280 for (i = 0; i < 3; i++)
2282 ir->deform[i][i] = dumdub[0][i];
2284 ir->deform[YY][XX] = dumdub[0][3];
2285 ir->deform[ZZ][XX] = dumdub[0][4];
2286 ir->deform[ZZ][YY] = dumdub[0][5];
2287 if (ir->epc != epcNO)
2289 for (i = 0; i < 3; i++)
2291 for (j = 0; j <= i; j++)
2293 if (ir->deform[i][j] != 0 && ir->compress[i][j] != 0)
2295 warning_error(wi, "A box element has deform set and compressibility > 0");
2299 for (i = 0; i < 3; i++)
2301 for (j = 0; j < i; j++)
2303 if (ir->deform[i][j] != 0)
2305 for (m = j; m < DIM; m++)
2307 if (ir->compress[m][j] != 0)
2309 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.");
2310 warning(wi, warn_buf);
2322 static int search_QMstring(const char *s, int ng, const char *gn[])
2324 /* same as normal search_string, but this one searches QM strings */
2327 for (i = 0; (i < ng); i++)
2329 if (gmx_strcasecmp(s, gn[i]) == 0)
2335 gmx_fatal(FARGS, "this QM method or basisset (%s) is not implemented\n!", s);
2339 } /* search_QMstring */
2341 /* We would like gn to be const as well, but C doesn't allow this */
2342 int search_string(const char *s, int ng, char *gn[])
2346 for (i = 0; (i < ng); i++)
2348 if (gmx_strcasecmp(s, gn[i]) == 0)
2355 "Group %s referenced in the .mdp file was not found in the index file.\n"
2356 "Group names must match either [moleculetype] names or custom index group\n"
2357 "names, in which case you must supply an index file to the '-n' option\n"
2364 static gmx_bool do_numbering(int natoms, gmx_groups_t *groups, int ng, char *ptrs[],
2365 t_blocka *block, char *gnames[],
2366 int gtype, int restnm,
2367 int grptp, gmx_bool bVerbose,
2370 unsigned short *cbuf;
2371 t_grps *grps = &(groups->grps[gtype]);
2372 int i, j, gid, aj, ognr, ntot = 0;
2375 char warn_buf[STRLEN];
2379 fprintf(debug, "Starting numbering %d groups of type %d\n", ng, gtype);
2382 title = gtypes[gtype];
2385 /* Mark all id's as not set */
2386 for (i = 0; (i < natoms); i++)
2391 snew(grps->nm_ind, ng+1); /* +1 for possible rest group */
2392 for (i = 0; (i < ng); i++)
2394 /* Lookup the group name in the block structure */
2395 gid = search_string(ptrs[i], block->nr, gnames);
2396 if ((grptp != egrptpONE) || (i == 0))
2398 grps->nm_ind[grps->nr++] = gid;
2402 fprintf(debug, "Found gid %d for group %s\n", gid, ptrs[i]);
2405 /* Now go over the atoms in the group */
2406 for (j = block->index[gid]; (j < block->index[gid+1]); j++)
2411 /* Range checking */
2412 if ((aj < 0) || (aj >= natoms))
2414 gmx_fatal(FARGS, "Invalid atom number %d in indexfile", aj);
2416 /* Lookup up the old group number */
2420 gmx_fatal(FARGS, "Atom %d in multiple %s groups (%d and %d)",
2421 aj+1, title, ognr+1, i+1);
2425 /* Store the group number in buffer */
2426 if (grptp == egrptpONE)
2439 /* Now check whether we have done all atoms */
2443 if (grptp == egrptpALL)
2445 gmx_fatal(FARGS, "%d atoms are not part of any of the %s groups",
2446 natoms-ntot, title);
2448 else if (grptp == egrptpPART)
2450 sprintf(warn_buf, "%d atoms are not part of any of the %s groups",
2451 natoms-ntot, title);
2452 warning_note(wi, warn_buf);
2454 /* Assign all atoms currently unassigned to a rest group */
2455 for (j = 0; (j < natoms); j++)
2457 if (cbuf[j] == NOGID)
2463 if (grptp != egrptpPART)
2468 "Making dummy/rest group for %s containing %d elements\n",
2469 title, natoms-ntot);
2471 /* Add group name "rest" */
2472 grps->nm_ind[grps->nr] = restnm;
2474 /* Assign the rest name to all atoms not currently assigned to a group */
2475 for (j = 0; (j < natoms); j++)
2477 if (cbuf[j] == NOGID)
2486 if (grps->nr == 1 && (ntot == 0 || ntot == natoms))
2488 /* All atoms are part of one (or no) group, no index required */
2489 groups->ngrpnr[gtype] = 0;
2490 groups->grpnr[gtype] = NULL;
2494 groups->ngrpnr[gtype] = natoms;
2495 snew(groups->grpnr[gtype], natoms);
2496 for (j = 0; (j < natoms); j++)
2498 groups->grpnr[gtype][j] = cbuf[j];
2504 return (bRest && grptp == egrptpPART);
2507 static void calc_nrdf(gmx_mtop_t *mtop, t_inputrec *ir, char **gnames)
2510 gmx_groups_t *groups;
2512 int natoms, ai, aj, i, j, d, g, imin, jmin;
2514 int *nrdf2, *na_vcm, na_tot;
2515 double *nrdf_tc, *nrdf_vcm, nrdf_uc, n_sub = 0;
2516 gmx_mtop_atomloop_all_t aloop;
2518 int mb, mol, ftype, as;
2519 gmx_molblock_t *molb;
2520 gmx_moltype_t *molt;
2523 * First calc 3xnr-atoms for each group
2524 * then subtract half a degree of freedom for each constraint
2526 * Only atoms and nuclei contribute to the degrees of freedom...
2531 groups = &mtop->groups;
2532 natoms = mtop->natoms;
2534 /* Allocate one more for a possible rest group */
2535 /* We need to sum degrees of freedom into doubles,
2536 * since floats give too low nrdf's above 3 million atoms.
2538 snew(nrdf_tc, groups->grps[egcTC].nr+1);
2539 snew(nrdf_vcm, groups->grps[egcVCM].nr+1);
2540 snew(na_vcm, groups->grps[egcVCM].nr+1);
2542 for (i = 0; i < groups->grps[egcTC].nr; i++)
2546 for (i = 0; i < groups->grps[egcVCM].nr+1; i++)
2551 snew(nrdf2, natoms);
2552 aloop = gmx_mtop_atomloop_all_init(mtop);
2553 while (gmx_mtop_atomloop_all_next(aloop, &i, &atom))
2556 if (atom->ptype == eptAtom || atom->ptype == eptNucleus)
2558 g = ggrpnr(groups, egcFREEZE, i);
2559 /* Double count nrdf for particle i */
2560 for (d = 0; d < DIM; d++)
2562 if (opts->nFreeze[g][d] == 0)
2567 nrdf_tc [ggrpnr(groups, egcTC, i)] += 0.5*nrdf2[i];
2568 nrdf_vcm[ggrpnr(groups, egcVCM, i)] += 0.5*nrdf2[i];
2573 for (mb = 0; mb < mtop->nmolblock; mb++)
2575 molb = &mtop->molblock[mb];
2576 molt = &mtop->moltype[molb->type];
2577 atom = molt->atoms.atom;
2578 for (mol = 0; mol < molb->nmol; mol++)
2580 for (ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++)
2582 ia = molt->ilist[ftype].iatoms;
2583 for (i = 0; i < molt->ilist[ftype].nr; )
2585 /* Subtract degrees of freedom for the constraints,
2586 * if the particles still have degrees of freedom left.
2587 * If one of the particles is a vsite or a shell, then all
2588 * constraint motion will go there, but since they do not
2589 * contribute to the constraints the degrees of freedom do not
2594 if (((atom[ia[1]].ptype == eptNucleus) ||
2595 (atom[ia[1]].ptype == eptAtom)) &&
2596 ((atom[ia[2]].ptype == eptNucleus) ||
2597 (atom[ia[2]].ptype == eptAtom)))
2615 imin = min(imin, nrdf2[ai]);
2616 jmin = min(jmin, nrdf2[aj]);
2619 nrdf_tc [ggrpnr(groups, egcTC, ai)] -= 0.5*imin;
2620 nrdf_tc [ggrpnr(groups, egcTC, aj)] -= 0.5*jmin;
2621 nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5*imin;
2622 nrdf_vcm[ggrpnr(groups, egcVCM, aj)] -= 0.5*jmin;
2624 ia += interaction_function[ftype].nratoms+1;
2625 i += interaction_function[ftype].nratoms+1;
2628 ia = molt->ilist[F_SETTLE].iatoms;
2629 for (i = 0; i < molt->ilist[F_SETTLE].nr; )
2631 /* Subtract 1 dof from every atom in the SETTLE */
2632 for (j = 0; j < 3; j++)
2635 imin = min(2, nrdf2[ai]);
2637 nrdf_tc [ggrpnr(groups, egcTC, ai)] -= 0.5*imin;
2638 nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5*imin;
2643 as += molt->atoms.nr;
2647 if (ir->ePull == epullCONSTRAINT)
2649 /* Correct nrdf for the COM constraints.
2650 * We correct using the TC and VCM group of the first atom
2651 * in the reference and pull group. If atoms in one pull group
2652 * belong to different TC or VCM groups it is anyhow difficult
2653 * to determine the optimal nrdf assignment.
2657 for (i = 0; i < pull->ncoord; i++)
2661 for (j = 0; j < 2; j++)
2663 const t_pull_group *pgrp;
2665 pgrp = &pull->group[pull->coord[i].group[j]];
2669 /* Subtract 1/2 dof from each group */
2671 nrdf_tc [ggrpnr(groups, egcTC, ai)] -= 0.5*imin;
2672 nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5*imin;
2673 if (nrdf_tc[ggrpnr(groups, egcTC, ai)] < 0)
2675 gmx_fatal(FARGS, "Center of mass pulling constraints caused the number of degrees of freedom for temperature coupling group %s to be negative", gnames[groups->grps[egcTC].nm_ind[ggrpnr(groups, egcTC, ai)]]);
2680 /* We need to subtract the whole DOF from group j=1 */
2687 if (ir->nstcomm != 0)
2689 /* Subtract 3 from the number of degrees of freedom in each vcm group
2690 * when com translation is removed and 6 when rotation is removed
2693 switch (ir->comm_mode)
2696 n_sub = ndof_com(ir);
2703 gmx_incons("Checking comm_mode");
2706 for (i = 0; i < groups->grps[egcTC].nr; i++)
2708 /* Count the number of atoms of TC group i for every VCM group */
2709 for (j = 0; j < groups->grps[egcVCM].nr+1; j++)
2714 for (ai = 0; ai < natoms; ai++)
2716 if (ggrpnr(groups, egcTC, ai) == i)
2718 na_vcm[ggrpnr(groups, egcVCM, ai)]++;
2722 /* Correct for VCM removal according to the fraction of each VCM
2723 * group present in this TC group.
2725 nrdf_uc = nrdf_tc[i];
2728 fprintf(debug, "T-group[%d] nrdf_uc = %g, n_sub = %g\n",
2732 for (j = 0; j < groups->grps[egcVCM].nr+1; j++)
2734 if (nrdf_vcm[j] > n_sub)
2736 nrdf_tc[i] += nrdf_uc*((double)na_vcm[j]/(double)na_tot)*
2737 (nrdf_vcm[j] - n_sub)/nrdf_vcm[j];
2741 fprintf(debug, " nrdf_vcm[%d] = %g, nrdf = %g\n",
2742 j, nrdf_vcm[j], nrdf_tc[i]);
2747 for (i = 0; (i < groups->grps[egcTC].nr); i++)
2749 opts->nrdf[i] = nrdf_tc[i];
2750 if (opts->nrdf[i] < 0)
2755 "Number of degrees of freedom in T-Coupling group %s is %.2f\n",
2756 gnames[groups->grps[egcTC].nm_ind[i]], opts->nrdf[i]);
2765 static void decode_cos(char *s, t_cosines *cosine)
2768 char format[STRLEN], f1[STRLEN];
2780 sscanf(t, "%d", &(cosine->n));
2787 snew(cosine->a, cosine->n);
2788 snew(cosine->phi, cosine->n);
2790 sprintf(format, "%%*d");
2791 for (i = 0; (i < cosine->n); i++)
2794 strcat(f1, "%lf%lf");
2795 if (sscanf(t, f1, &a, &phi) < 2)
2797 gmx_fatal(FARGS, "Invalid input for electric field shift: '%s'", t);
2800 cosine->phi[i] = phi;
2801 strcat(format, "%*lf%*lf");
2808 static gmx_bool do_egp_flag(t_inputrec *ir, gmx_groups_t *groups,
2809 const char *option, const char *val, int flag)
2811 /* The maximum number of energy group pairs would be MAXPTR*(MAXPTR+1)/2.
2812 * But since this is much larger than STRLEN, such a line can not be parsed.
2813 * The real maximum is the number of names that fit in a string: STRLEN/2.
2815 #define EGP_MAX (STRLEN/2)
2816 int nelem, i, j, k, nr;
2817 char *names[EGP_MAX];
2821 gnames = groups->grpname;
2823 nelem = str_nelem(val, EGP_MAX, names);
2826 gmx_fatal(FARGS, "The number of groups for %s is odd", option);
2828 nr = groups->grps[egcENER].nr;
2830 for (i = 0; i < nelem/2; i++)
2834 gmx_strcasecmp(names[2*i], *(gnames[groups->grps[egcENER].nm_ind[j]])))
2840 gmx_fatal(FARGS, "%s in %s is not an energy group\n",
2841 names[2*i], option);
2845 gmx_strcasecmp(names[2*i+1], *(gnames[groups->grps[egcENER].nm_ind[k]])))
2851 gmx_fatal(FARGS, "%s in %s is not an energy group\n",
2852 names[2*i+1], option);
2854 if ((j < nr) && (k < nr))
2856 ir->opts.egp_flags[nr*j+k] |= flag;
2857 ir->opts.egp_flags[nr*k+j] |= flag;
2865 void do_index(const char* mdparin, const char *ndx,
2868 t_inputrec *ir, rvec *v,
2872 gmx_groups_t *groups;
2876 char warnbuf[STRLEN], **gnames;
2877 int nr, ntcg, ntau_t, nref_t, nacc, nofg, nSA, nSA_points, nSA_time, nSA_temp;
2880 int nacg, nfreeze, nfrdim, nenergy, nvcm, nuser;
2881 char *ptr1[MAXPTR], *ptr2[MAXPTR], *ptr3[MAXPTR];
2882 int i, j, k, restnm;
2884 gmx_bool bExcl, bTable, bSetTCpar, bAnneal, bRest;
2885 int nQMmethod, nQMbasis, nQMcharge, nQMmult, nbSH, nCASorb, nCASelec,
2886 nSAon, nSAoff, nSAsteps, nQMg, nbOPT, nbTS;
2887 char warn_buf[STRLEN];
2891 fprintf(stderr, "processing index file...\n");
2897 snew(grps->index, 1);
2899 atoms_all = gmx_mtop_global_atoms(mtop);
2900 analyse(&atoms_all, grps, &gnames, FALSE, TRUE);
2901 free_t_atoms(&atoms_all, FALSE);
2905 grps = init_index(ndx, &gnames);
2908 groups = &mtop->groups;
2909 natoms = mtop->natoms;
2910 symtab = &mtop->symtab;
2912 snew(groups->grpname, grps->nr+1);
2914 for (i = 0; (i < grps->nr); i++)
2916 groups->grpname[i] = put_symtab(symtab, gnames[i]);
2918 groups->grpname[i] = put_symtab(symtab, "rest");
2920 srenew(gnames, grps->nr+1);
2921 gnames[restnm] = *(groups->grpname[i]);
2922 groups->ngrpname = grps->nr+1;
2924 set_warning_line(wi, mdparin, -1);
2926 ntau_t = str_nelem(tau_t, MAXPTR, ptr1);
2927 nref_t = str_nelem(ref_t, MAXPTR, ptr2);
2928 ntcg = str_nelem(tcgrps, MAXPTR, ptr3);
2929 if ((ntau_t != ntcg) || (nref_t != ntcg))
2931 gmx_fatal(FARGS, "Invalid T coupling input: %d groups, %d ref-t values and "
2932 "%d tau-t values", ntcg, nref_t, ntau_t);
2935 bSetTCpar = (ir->etc || EI_SD(ir->eI) || ir->eI == eiBD || EI_TPI(ir->eI));
2936 do_numbering(natoms, groups, ntcg, ptr3, grps, gnames, egcTC,
2937 restnm, bSetTCpar ? egrptpALL : egrptpALL_GENREST, bVerbose, wi);
2938 nr = groups->grps[egcTC].nr;
2940 snew(ir->opts.nrdf, nr);
2941 snew(ir->opts.tau_t, nr);
2942 snew(ir->opts.ref_t, nr);
2943 if (ir->eI == eiBD && ir->bd_fric == 0)
2945 fprintf(stderr, "bd-fric=0, so tau-t will be used as the inverse friction constant(s)\n");
2952 gmx_fatal(FARGS, "Not enough ref-t and tau-t values!");
2956 for (i = 0; (i < nr); i++)
2958 ir->opts.tau_t[i] = strtod(ptr1[i], NULL);
2959 if ((ir->eI == eiBD || ir->eI == eiSD2) && ir->opts.tau_t[i] <= 0)
2961 sprintf(warn_buf, "With integrator %s tau-t should be larger than 0", ei_names[ir->eI]);
2962 warning_error(wi, warn_buf);
2965 if (ir->etc != etcVRESCALE && ir->opts.tau_t[i] == 0)
2967 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.");
2970 if (ir->opts.tau_t[i] >= 0)
2972 tau_min = min(tau_min, ir->opts.tau_t[i]);
2975 if (ir->etc != etcNO && ir->nsttcouple == -1)
2977 ir->nsttcouple = ir_optimal_nsttcouple(ir);
2982 if ((ir->etc == etcNOSEHOOVER) && (ir->epc == epcBERENDSEN))
2984 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");
2986 if ((ir->epc == epcMTTK) && (ir->etc > etcNO))
2988 if (ir->nstpcouple != ir->nsttcouple)
2990 int mincouple = min(ir->nstpcouple, ir->nsttcouple);
2991 ir->nstpcouple = ir->nsttcouple = mincouple;
2992 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);
2993 warning_note(wi, warn_buf);
2997 /* velocity verlet with averaged kinetic energy KE = 0.5*(v(t+1/2) - v(t-1/2)) is implemented
2998 primarily for testing purposes, and does not work with temperature coupling other than 1 */
3000 if (ETC_ANDERSEN(ir->etc))
3002 if (ir->nsttcouple != 1)
3005 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");
3006 warning_note(wi, warn_buf);
3009 nstcmin = tcouple_min_integration_steps(ir->etc);
3012 if (tau_min/(ir->delta_t*ir->nsttcouple) < nstcmin)
3014 sprintf(warn_buf, "For proper integration of the %s thermostat, tau-t (%g) should be at least %d times larger than nsttcouple*dt (%g)",
3015 ETCOUPLTYPE(ir->etc),
3017 ir->nsttcouple*ir->delta_t);
3018 warning(wi, warn_buf);
3021 for (i = 0; (i < nr); i++)
3023 ir->opts.ref_t[i] = strtod(ptr2[i], NULL);
3024 if (ir->opts.ref_t[i] < 0)
3026 gmx_fatal(FARGS, "ref-t for group %d negative", i);
3029 /* set the lambda mc temperature to the md integrator temperature (which should be defined
3030 if we are in this conditional) if mc_temp is negative */
3031 if (ir->expandedvals->mc_temp < 0)
3033 ir->expandedvals->mc_temp = ir->opts.ref_t[0]; /*for now, set to the first reft */
3037 /* Simulated annealing for each group. There are nr groups */
3038 nSA = str_nelem(anneal, MAXPTR, ptr1);
3039 if (nSA == 1 && (ptr1[0][0] == 'n' || ptr1[0][0] == 'N'))
3043 if (nSA > 0 && nSA != nr)
3045 gmx_fatal(FARGS, "Not enough annealing values: %d (for %d groups)\n", nSA, nr);
3049 snew(ir->opts.annealing, nr);
3050 snew(ir->opts.anneal_npoints, nr);
3051 snew(ir->opts.anneal_time, nr);
3052 snew(ir->opts.anneal_temp, nr);
3053 for (i = 0; i < nr; i++)
3055 ir->opts.annealing[i] = eannNO;
3056 ir->opts.anneal_npoints[i] = 0;
3057 ir->opts.anneal_time[i] = NULL;
3058 ir->opts.anneal_temp[i] = NULL;
3063 for (i = 0; i < nr; i++)
3065 if (ptr1[i][0] == 'n' || ptr1[i][0] == 'N')
3067 ir->opts.annealing[i] = eannNO;
3069 else if (ptr1[i][0] == 's' || ptr1[i][0] == 'S')
3071 ir->opts.annealing[i] = eannSINGLE;
3074 else if (ptr1[i][0] == 'p' || ptr1[i][0] == 'P')
3076 ir->opts.annealing[i] = eannPERIODIC;
3082 /* Read the other fields too */
3083 nSA_points = str_nelem(anneal_npoints, MAXPTR, ptr1);
3084 if (nSA_points != nSA)
3086 gmx_fatal(FARGS, "Found %d annealing-npoints values for %d groups\n", nSA_points, nSA);
3088 for (k = 0, i = 0; i < nr; i++)
3090 ir->opts.anneal_npoints[i] = strtol(ptr1[i], NULL, 10);
3091 if (ir->opts.anneal_npoints[i] == 1)
3093 gmx_fatal(FARGS, "Please specify at least a start and an end point for annealing\n");
3095 snew(ir->opts.anneal_time[i], ir->opts.anneal_npoints[i]);
3096 snew(ir->opts.anneal_temp[i], ir->opts.anneal_npoints[i]);
3097 k += ir->opts.anneal_npoints[i];
3100 nSA_time = str_nelem(anneal_time, MAXPTR, ptr1);
3103 gmx_fatal(FARGS, "Found %d annealing-time values, wanter %d\n", nSA_time, k);
3105 nSA_temp = str_nelem(anneal_temp, MAXPTR, ptr2);
3108 gmx_fatal(FARGS, "Found %d annealing-temp values, wanted %d\n", nSA_temp, k);
3111 for (i = 0, k = 0; i < nr; i++)
3114 for (j = 0; j < ir->opts.anneal_npoints[i]; j++)
3116 ir->opts.anneal_time[i][j] = strtod(ptr1[k], NULL);
3117 ir->opts.anneal_temp[i][j] = strtod(ptr2[k], NULL);
3120 if (ir->opts.anneal_time[i][0] > (ir->init_t+GMX_REAL_EPS))
3122 gmx_fatal(FARGS, "First time point for annealing > init_t.\n");
3128 if (ir->opts.anneal_time[i][j] < ir->opts.anneal_time[i][j-1])
3130 gmx_fatal(FARGS, "Annealing timepoints out of order: t=%f comes after t=%f\n",
3131 ir->opts.anneal_time[i][j], ir->opts.anneal_time[i][j-1]);
3134 if (ir->opts.anneal_temp[i][j] < 0)
3136 gmx_fatal(FARGS, "Found negative temperature in annealing: %f\n", ir->opts.anneal_temp[i][j]);
3141 /* Print out some summary information, to make sure we got it right */
3142 for (i = 0, k = 0; i < nr; i++)
3144 if (ir->opts.annealing[i] != eannNO)
3146 j = groups->grps[egcTC].nm_ind[i];
3147 fprintf(stderr, "Simulated annealing for group %s: %s, %d timepoints\n",
3148 *(groups->grpname[j]), eann_names[ir->opts.annealing[i]],
3149 ir->opts.anneal_npoints[i]);
3150 fprintf(stderr, "Time (ps) Temperature (K)\n");
3151 /* All terms except the last one */
3152 for (j = 0; j < (ir->opts.anneal_npoints[i]-1); j++)
3154 fprintf(stderr, "%9.1f %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3157 /* Finally the last one */
3158 j = ir->opts.anneal_npoints[i]-1;
3159 if (ir->opts.annealing[i] == eannSINGLE)
3161 fprintf(stderr, "%9.1f- %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3165 fprintf(stderr, "%9.1f %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3166 if (fabs(ir->opts.anneal_temp[i][j]-ir->opts.anneal_temp[i][0]) > GMX_REAL_EPS)
3168 warning_note(wi, "There is a temperature jump when your annealing loops back.\n");
3177 if (ir->ePull != epullNO)
3179 make_pull_groups(ir->pull, pull_grp, grps, gnames);
3181 make_pull_coords(ir->pull);
3186 make_rotation_groups(ir->rot, rot_grp, grps, gnames);
3189 nacc = str_nelem(acc, MAXPTR, ptr1);
3190 nacg = str_nelem(accgrps, MAXPTR, ptr2);
3191 if (nacg*DIM != nacc)
3193 gmx_fatal(FARGS, "Invalid Acceleration input: %d groups and %d acc. values",
3196 do_numbering(natoms, groups, nacg, ptr2, grps, gnames, egcACC,
3197 restnm, egrptpALL_GENREST, bVerbose, wi);
3198 nr = groups->grps[egcACC].nr;
3199 snew(ir->opts.acc, nr);
3200 ir->opts.ngacc = nr;
3202 for (i = k = 0; (i < nacg); i++)
3204 for (j = 0; (j < DIM); j++, k++)
3206 ir->opts.acc[i][j] = strtod(ptr1[k], NULL);
3209 for (; (i < nr); i++)
3211 for (j = 0; (j < DIM); j++)
3213 ir->opts.acc[i][j] = 0;
3217 nfrdim = str_nelem(frdim, MAXPTR, ptr1);
3218 nfreeze = str_nelem(freeze, MAXPTR, ptr2);
3219 if (nfrdim != DIM*nfreeze)
3221 gmx_fatal(FARGS, "Invalid Freezing input: %d groups and %d freeze values",
3224 do_numbering(natoms, groups, nfreeze, ptr2, grps, gnames, egcFREEZE,
3225 restnm, egrptpALL_GENREST, bVerbose, wi);
3226 nr = groups->grps[egcFREEZE].nr;
3227 ir->opts.ngfrz = nr;
3228 snew(ir->opts.nFreeze, nr);
3229 for (i = k = 0; (i < nfreeze); i++)
3231 for (j = 0; (j < DIM); j++, k++)
3233 ir->opts.nFreeze[i][j] = (gmx_strncasecmp(ptr1[k], "Y", 1) == 0);
3234 if (!ir->opts.nFreeze[i][j])
3236 if (gmx_strncasecmp(ptr1[k], "N", 1) != 0)
3238 sprintf(warnbuf, "Please use Y(ES) or N(O) for freezedim only "
3239 "(not %s)", ptr1[k]);
3240 warning(wi, warn_buf);
3245 for (; (i < nr); i++)
3247 for (j = 0; (j < DIM); j++)
3249 ir->opts.nFreeze[i][j] = 0;
3253 nenergy = str_nelem(energy, MAXPTR, ptr1);
3254 do_numbering(natoms, groups, nenergy, ptr1, grps, gnames, egcENER,
3255 restnm, egrptpALL_GENREST, bVerbose, wi);
3256 add_wall_energrps(groups, ir->nwall, symtab);
3257 ir->opts.ngener = groups->grps[egcENER].nr;
3258 nvcm = str_nelem(vcm, MAXPTR, ptr1);
3260 do_numbering(natoms, groups, nvcm, ptr1, grps, gnames, egcVCM,
3261 restnm, nvcm == 0 ? egrptpALL_GENREST : egrptpPART, bVerbose, wi);
3264 warning(wi, "Some atoms are not part of any center of mass motion removal group.\n"
3265 "This may lead to artifacts.\n"
3266 "In most cases one should use one group for the whole system.");
3269 /* Now we have filled the freeze struct, so we can calculate NRDF */
3270 calc_nrdf(mtop, ir, gnames);
3276 /* Must check per group! */
3277 for (i = 0; (i < ir->opts.ngtc); i++)
3279 ntot += ir->opts.nrdf[i];
3281 if (ntot != (DIM*natoms))
3283 fac = sqrt(ntot/(DIM*natoms));
3286 fprintf(stderr, "Scaling velocities by a factor of %.3f to account for constraints\n"
3287 "and removal of center of mass motion\n", fac);
3289 for (i = 0; (i < natoms); i++)
3291 svmul(fac, v[i], v[i]);
3296 nuser = str_nelem(user1, MAXPTR, ptr1);
3297 do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcUser1,
3298 restnm, egrptpALL_GENREST, bVerbose, wi);
3299 nuser = str_nelem(user2, MAXPTR, ptr1);
3300 do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcUser2,
3301 restnm, egrptpALL_GENREST, bVerbose, wi);
3302 nuser = str_nelem(xtc_grps, MAXPTR, ptr1);
3303 do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcXTC,
3304 restnm, egrptpONE, bVerbose, wi);
3305 nofg = str_nelem(orirefitgrp, MAXPTR, ptr1);
3306 do_numbering(natoms, groups, nofg, ptr1, grps, gnames, egcORFIT,
3307 restnm, egrptpALL_GENREST, bVerbose, wi);
3309 /* QMMM input processing */
3310 nQMg = str_nelem(QMMM, MAXPTR, ptr1);
3311 nQMmethod = str_nelem(QMmethod, MAXPTR, ptr2);
3312 nQMbasis = str_nelem(QMbasis, MAXPTR, ptr3);
3313 if ((nQMmethod != nQMg) || (nQMbasis != nQMg))
3315 gmx_fatal(FARGS, "Invalid QMMM input: %d groups %d basissets"
3316 " and %d methods\n", nQMg, nQMbasis, nQMmethod);
3318 /* group rest, if any, is always MM! */
3319 do_numbering(natoms, groups, nQMg, ptr1, grps, gnames, egcQMMM,
3320 restnm, egrptpALL_GENREST, bVerbose, wi);
3321 nr = nQMg; /*atoms->grps[egcQMMM].nr;*/
3322 ir->opts.ngQM = nQMg;
3323 snew(ir->opts.QMmethod, nr);
3324 snew(ir->opts.QMbasis, nr);
3325 for (i = 0; i < nr; i++)
3327 /* input consists of strings: RHF CASSCF PM3 .. These need to be
3328 * converted to the corresponding enum in names.c
3330 ir->opts.QMmethod[i] = search_QMstring(ptr2[i], eQMmethodNR,
3332 ir->opts.QMbasis[i] = search_QMstring(ptr3[i], eQMbasisNR,
3336 nQMmult = str_nelem(QMmult, MAXPTR, ptr1);
3337 nQMcharge = str_nelem(QMcharge, MAXPTR, ptr2);
3338 nbSH = str_nelem(bSH, MAXPTR, ptr3);
3339 snew(ir->opts.QMmult, nr);
3340 snew(ir->opts.QMcharge, nr);
3341 snew(ir->opts.bSH, nr);
3343 for (i = 0; i < nr; i++)
3345 ir->opts.QMmult[i] = strtol(ptr1[i], NULL, 10);
3346 ir->opts.QMcharge[i] = strtol(ptr2[i], NULL, 10);
3347 ir->opts.bSH[i] = (gmx_strncasecmp(ptr3[i], "Y", 1) == 0);
3350 nCASelec = str_nelem(CASelectrons, MAXPTR, ptr1);
3351 nCASorb = str_nelem(CASorbitals, MAXPTR, ptr2);
3352 snew(ir->opts.CASelectrons, nr);
3353 snew(ir->opts.CASorbitals, nr);
3354 for (i = 0; i < nr; i++)
3356 ir->opts.CASelectrons[i] = strtol(ptr1[i], NULL, 10);
3357 ir->opts.CASorbitals[i] = strtol(ptr2[i], NULL, 10);
3359 /* special optimization options */
3361 nbOPT = str_nelem(bOPT, MAXPTR, ptr1);
3362 nbTS = str_nelem(bTS, MAXPTR, ptr2);
3363 snew(ir->opts.bOPT, nr);
3364 snew(ir->opts.bTS, nr);
3365 for (i = 0; i < nr; i++)
3367 ir->opts.bOPT[i] = (gmx_strncasecmp(ptr1[i], "Y", 1) == 0);
3368 ir->opts.bTS[i] = (gmx_strncasecmp(ptr2[i], "Y", 1) == 0);
3370 nSAon = str_nelem(SAon, MAXPTR, ptr1);
3371 nSAoff = str_nelem(SAoff, MAXPTR, ptr2);
3372 nSAsteps = str_nelem(SAsteps, MAXPTR, ptr3);
3373 snew(ir->opts.SAon, nr);
3374 snew(ir->opts.SAoff, nr);
3375 snew(ir->opts.SAsteps, nr);
3377 for (i = 0; i < nr; i++)
3379 ir->opts.SAon[i] = strtod(ptr1[i], NULL);
3380 ir->opts.SAoff[i] = strtod(ptr2[i], NULL);
3381 ir->opts.SAsteps[i] = strtol(ptr3[i], NULL, 10);
3383 /* end of QMMM input */
3387 for (i = 0; (i < egcNR); i++)
3389 fprintf(stderr, "%-16s has %d element(s):", gtypes[i], groups->grps[i].nr);
3390 for (j = 0; (j < groups->grps[i].nr); j++)
3392 fprintf(stderr, " %s", *(groups->grpname[groups->grps[i].nm_ind[j]]));
3394 fprintf(stderr, "\n");
3398 nr = groups->grps[egcENER].nr;
3399 snew(ir->opts.egp_flags, nr*nr);
3401 bExcl = do_egp_flag(ir, groups, "energygrp-excl", egpexcl, EGP_EXCL);
3402 if (bExcl && ir->cutoff_scheme == ecutsVERLET)
3404 warning_error(wi, "Energy group exclusions are not (yet) implemented for the Verlet scheme");
3406 if (bExcl && EEL_FULL(ir->coulombtype))
3408 warning(wi, "Can not exclude the lattice Coulomb energy between energy groups");
3411 bTable = do_egp_flag(ir, groups, "energygrp-table", egptable, EGP_TABLE);
3412 if (bTable && !(ir->vdwtype == evdwUSER) &&
3413 !(ir->coulombtype == eelUSER) && !(ir->coulombtype == eelPMEUSER) &&
3414 !(ir->coulombtype == eelPMEUSERSWITCH))
3416 gmx_fatal(FARGS, "Can only have energy group pair tables in combination with user tables for VdW and/or Coulomb");
3419 decode_cos(efield_x, &(ir->ex[XX]));
3420 decode_cos(efield_xt, &(ir->et[XX]));
3421 decode_cos(efield_y, &(ir->ex[YY]));
3422 decode_cos(efield_yt, &(ir->et[YY]));
3423 decode_cos(efield_z, &(ir->ex[ZZ]));
3424 decode_cos(efield_zt, &(ir->et[ZZ]));
3428 do_adress_index(ir->adress, groups, gnames, &(ir->opts), wi);
3431 for (i = 0; (i < grps->nr); i++)
3443 static void check_disre(gmx_mtop_t *mtop)
3445 gmx_ffparams_t *ffparams;
3446 t_functype *functype;
3448 int i, ndouble, ftype;
3449 int label, old_label;
3451 if (gmx_mtop_ftype_count(mtop, F_DISRES) > 0)
3453 ffparams = &mtop->ffparams;
3454 functype = ffparams->functype;
3455 ip = ffparams->iparams;
3458 for (i = 0; i < ffparams->ntypes; i++)
3460 ftype = functype[i];
3461 if (ftype == F_DISRES)
3463 label = ip[i].disres.label;
3464 if (label == old_label)
3466 fprintf(stderr, "Distance restraint index %d occurs twice\n", label);
3474 gmx_fatal(FARGS, "Found %d double distance restraint indices,\n"
3475 "probably the parameters for multiple pairs in one restraint "
3476 "are not identical\n", ndouble);
3481 static gmx_bool absolute_reference(t_inputrec *ir, gmx_mtop_t *sys,
3482 gmx_bool posres_only,
3486 gmx_mtop_ilistloop_t iloop;
3496 for (d = 0; d < DIM; d++)
3498 AbsRef[d] = (d < ndof_com(ir) ? 0 : 1);
3500 /* Check for freeze groups */
3501 for (g = 0; g < ir->opts.ngfrz; g++)
3503 for (d = 0; d < DIM; d++)
3505 if (ir->opts.nFreeze[g][d] != 0)
3513 /* Check for position restraints */
3514 iloop = gmx_mtop_ilistloop_init(sys);
3515 while (gmx_mtop_ilistloop_next(iloop, &ilist, &nmol))
3518 (AbsRef[XX] == 0 || AbsRef[YY] == 0 || AbsRef[ZZ] == 0))
3520 for (i = 0; i < ilist[F_POSRES].nr; i += 2)
3522 pr = &sys->ffparams.iparams[ilist[F_POSRES].iatoms[i]];
3523 for (d = 0; d < DIM; d++)
3525 if (pr->posres.fcA[d] != 0)
3531 for (i = 0; i < ilist[F_FBPOSRES].nr; i += 2)
3533 /* Check for flat-bottom posres */
3534 pr = &sys->ffparams.iparams[ilist[F_FBPOSRES].iatoms[i]];
3535 if (pr->fbposres.k != 0)
3537 switch (pr->fbposres.geom)
3539 case efbposresSPHERE:
3540 AbsRef[XX] = AbsRef[YY] = AbsRef[ZZ] = 1;
3542 case efbposresCYLINDER:
3543 AbsRef[XX] = AbsRef[YY] = 1;
3545 case efbposresX: /* d=XX */
3546 case efbposresY: /* d=YY */
3547 case efbposresZ: /* d=ZZ */
3548 d = pr->fbposres.geom - efbposresX;
3552 gmx_fatal(FARGS, " Invalid geometry for flat-bottom position restraint.\n"
3553 "Expected nr between 1 and %d. Found %d\n", efbposresNR-1,
3561 return (AbsRef[XX] != 0 && AbsRef[YY] != 0 && AbsRef[ZZ] != 0);
3565 check_combination_rule_differences(const gmx_mtop_t *mtop, int state,
3566 gmx_bool *bC6ParametersWorkWithGeometricRules,
3567 gmx_bool *bC6ParametersWorkWithLBRules,
3568 gmx_bool *bLBRulesPossible)
3570 int ntypes, tpi, tpj, thisLBdiff, thisgeomdiff;
3573 #if (defined SIZEOF_LONG_LONG_INT) && (SIZEOF_LONG_LONG_INT >= 8)
3574 long long int npair, npair_ij;
3576 double npair, npair_ij;
3578 double geometricdiff, LBdiff;
3579 double c6i, c6j, c12i, c12j;
3580 double c6, c6_geometric, c6_LB;
3581 double sigmai, sigmaj, epsi, epsj;
3582 gmx_bool bCanDoLBRules, bCanDoGeometricRules;
3585 /* A tolerance of 1e-5 seems reasonable for (possibly hand-typed)
3586 * force-field floating point parameters.
3589 ptr = getenv("GMX_LJCOMB_TOL");
3594 sscanf(ptr, "%lf", &dbl);
3598 *bC6ParametersWorkWithLBRules = TRUE;
3599 *bC6ParametersWorkWithGeometricRules = TRUE;
3600 bCanDoLBRules = TRUE;
3601 bCanDoGeometricRules = TRUE;
3603 ntypes = mtop->ffparams.atnr;
3604 snew(typecount, ntypes);
3605 gmx_mtop_count_atomtypes(mtop, state, typecount);
3606 geometricdiff = LBdiff = 0.0;
3607 *bLBRulesPossible = TRUE;
3608 for (tpi = 0; tpi < ntypes; ++tpi)
3610 c6i = mtop->ffparams.iparams[(ntypes + 1) * tpi].lj.c6;
3611 c12i = mtop->ffparams.iparams[(ntypes + 1) * tpi].lj.c12;
3612 for (tpj = tpi; tpj < ntypes; ++tpj)
3614 c6j = mtop->ffparams.iparams[(ntypes + 1) * tpj].lj.c6;
3615 c12j = mtop->ffparams.iparams[(ntypes + 1) * tpj].lj.c12;
3616 c6 = mtop->ffparams.iparams[ntypes * tpi + tpj].lj.c6;
3617 c6_geometric = sqrt(c6i * c6j);
3618 if (!gmx_numzero(c6_geometric))
3620 if (!gmx_numzero(c12i) && !gmx_numzero(c12j))
3622 sigmai = pow(c12i / c6i, 1.0/6.0);
3623 sigmaj = pow(c12j / c6j, 1.0/6.0);
3624 epsi = c6i * c6i /(4.0 * c12i);
3625 epsj = c6j * c6j /(4.0 * c12j);
3626 c6_LB = 4.0 * pow(epsi * epsj, 1.0/2.0) * pow(0.5 * (sigmai + sigmaj), 6);
3630 *bLBRulesPossible = FALSE;
3631 c6_LB = c6_geometric;
3633 bCanDoLBRules = gmx_within_tol(c6_LB, c6, tol);
3636 if (FALSE == bCanDoLBRules)
3638 *bC6ParametersWorkWithLBRules = FALSE;
3641 bCanDoGeometricRules = gmx_within_tol(c6_geometric, c6, tol);
3643 if (FALSE == bCanDoGeometricRules)
3645 *bC6ParametersWorkWithGeometricRules = FALSE;
3653 check_combination_rules(const t_inputrec *ir, const gmx_mtop_t *mtop,
3657 gmx_bool bLBRulesPossible, bC6ParametersWorkWithGeometricRules, bC6ParametersWorkWithLBRules;
3659 check_combination_rule_differences(mtop, 0,
3660 &bC6ParametersWorkWithGeometricRules,
3661 &bC6ParametersWorkWithLBRules,
3663 if (ir->ljpme_combination_rule == eljpmeLB)
3665 if (FALSE == bC6ParametersWorkWithLBRules || FALSE == bLBRulesPossible)
3667 warning(wi, "You are using arithmetic-geometric combination rules "
3668 "in LJ-PME, but your non-bonded C6 parameters do not "
3669 "follow these rules.");
3674 if (FALSE == bC6ParametersWorkWithGeometricRules)
3676 if (ir->eDispCorr != edispcNO)
3678 warning_note(wi, "You are using geometric combination rules in "
3679 "LJ-PME, but your non-bonded C6 parameters do "
3680 "not follow these rules. "
3681 "If your force field uses Lorentz-Berthelot combination rules this "
3682 "will introduce small errors in the forces and energies in "
3683 "your simulations. Dispersion correction will correct total "
3684 "energy and/or pressure, but not forces or surface tensions. "
3685 "Please check the LJ-PME section in the manual "
3686 "before proceeding further.");
3690 warning_note(wi, "You are using geometric combination rules in "
3691 "LJ-PME, but your non-bonded C6 parameters do "
3692 "not follow these rules. "
3693 "If your force field uses Lorentz-Berthelot combination rules this "
3694 "will introduce small errors in the forces and energies in "
3695 "your simulations. Consider using dispersion correction "
3696 "for the total energy and pressure. "
3697 "Please check the LJ-PME section in the manual "
3698 "before proceeding further.");
3704 void triple_check(const char *mdparin, t_inputrec *ir, gmx_mtop_t *sys,
3708 int i, m, c, nmol, npct;
3709 gmx_bool bCharge, bAcc;
3710 real gdt_max, *mgrp, mt;
3712 gmx_mtop_atomloop_block_t aloopb;
3713 gmx_mtop_atomloop_all_t aloop;
3716 char warn_buf[STRLEN];
3718 set_warning_line(wi, mdparin, -1);
3720 if (EI_DYNAMICS(ir->eI) && !EI_SD(ir->eI) && ir->eI != eiBD &&
3721 ir->comm_mode == ecmNO &&
3722 !(absolute_reference(ir, sys, FALSE, AbsRef) || ir->nsteps <= 10))
3724 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");
3727 /* Check for pressure coupling with absolute position restraints */
3728 if (ir->epc != epcNO && ir->refcoord_scaling == erscNO)
3730 absolute_reference(ir, sys, TRUE, AbsRef);
3732 for (m = 0; m < DIM; m++)
3734 if (AbsRef[m] && norm2(ir->compress[m]) > 0)
3736 warning(wi, "You are using pressure coupling with absolute position restraints, this will give artifacts. Use the refcoord_scaling option.");
3744 aloopb = gmx_mtop_atomloop_block_init(sys);
3745 while (gmx_mtop_atomloop_block_next(aloopb, &atom, &nmol))
3747 if (atom->q != 0 || atom->qB != 0)
3755 if (EEL_FULL(ir->coulombtype))
3758 "You are using full electrostatics treatment %s for a system without charges.\n"
3759 "This costs a lot of performance for just processing zeros, consider using %s instead.\n",
3760 EELTYPE(ir->coulombtype), EELTYPE(eelCUT));
3761 warning(wi, err_buf);
3766 if (ir->coulombtype == eelCUT && ir->rcoulomb > 0 && !ir->implicit_solvent)
3769 "You are using a plain Coulomb cut-off, which might produce artifacts.\n"
3770 "You might want to consider using %s electrostatics.\n",
3772 warning_note(wi, err_buf);
3776 /* Check if combination rules used in LJ-PME are the same as in the force field */
3777 if (EVDW_PME(ir->vdwtype))
3779 check_combination_rules(ir, sys, wi);
3782 /* Generalized reaction field */
3783 if (ir->opts.ngtc == 0)
3785 sprintf(err_buf, "No temperature coupling while using coulombtype %s",
3787 CHECK(ir->coulombtype == eelGRF);
3791 sprintf(err_buf, "When using coulombtype = %s"
3792 " ref-t for temperature coupling should be > 0",
3794 CHECK((ir->coulombtype == eelGRF) && (ir->opts.ref_t[0] <= 0));
3797 if (ir->eI == eiSD1 &&
3798 (gmx_mtop_ftype_count(sys, F_CONSTR) > 0 ||
3799 gmx_mtop_ftype_count(sys, F_SETTLE) > 0))
3801 sprintf(warn_buf, "With constraints integrator %s is less accurate, consider using %s instead", ei_names[ir->eI], ei_names[eiSD2]);
3802 warning_note(wi, warn_buf);
3806 for (i = 0; (i < sys->groups.grps[egcACC].nr); i++)
3808 for (m = 0; (m < DIM); m++)
3810 if (fabs(ir->opts.acc[i][m]) > 1e-6)
3819 snew(mgrp, sys->groups.grps[egcACC].nr);
3820 aloop = gmx_mtop_atomloop_all_init(sys);
3821 while (gmx_mtop_atomloop_all_next(aloop, &i, &atom))
3823 mgrp[ggrpnr(&sys->groups, egcACC, i)] += atom->m;
3826 for (i = 0; (i < sys->groups.grps[egcACC].nr); i++)
3828 for (m = 0; (m < DIM); m++)
3830 acc[m] += ir->opts.acc[i][m]*mgrp[i];
3834 for (m = 0; (m < DIM); m++)
3836 if (fabs(acc[m]) > 1e-6)
3838 const char *dim[DIM] = { "X", "Y", "Z" };
3840 "Net Acceleration in %s direction, will %s be corrected\n",
3841 dim[m], ir->nstcomm != 0 ? "" : "not");
3842 if (ir->nstcomm != 0 && m < ndof_com(ir))
3845 for (i = 0; (i < sys->groups.grps[egcACC].nr); i++)
3847 ir->opts.acc[i][m] -= acc[m];
3855 if (ir->efep != efepNO && ir->fepvals->sc_alpha != 0 &&
3856 !gmx_within_tol(sys->ffparams.reppow, 12.0, 10*GMX_DOUBLE_EPS))
3858 gmx_fatal(FARGS, "Soft-core interactions are only supported with VdW repulsion power 12");
3861 if (ir->ePull != epullNO)
3863 gmx_bool bPullAbsoluteRef;
3865 bPullAbsoluteRef = FALSE;
3866 for (i = 0; i < ir->pull->ncoord; i++)
3868 bPullAbsoluteRef = bPullAbsoluteRef ||
3869 ir->pull->coord[i].group[0] == 0 ||
3870 ir->pull->coord[i].group[1] == 0;
3872 if (bPullAbsoluteRef)
3874 absolute_reference(ir, sys, FALSE, AbsRef);
3875 for (m = 0; m < DIM; m++)
3877 if (ir->pull->dim[m] && !AbsRef[m])
3879 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.");
3885 if (ir->pull->eGeom == epullgDIRPBC)
3887 for (i = 0; i < 3; i++)
3889 for (m = 0; m <= i; m++)
3891 if ((ir->epc != epcNO && ir->compress[i][m] != 0) ||
3892 ir->deform[i][m] != 0)
3894 for (c = 0; c < ir->pull->ncoord; c++)
3896 if (ir->pull->coord[c].vec[m] != 0)
3898 gmx_fatal(FARGS, "Can not have dynamic box while using pull geometry '%s' (dim %c)", EPULLGEOM(ir->pull->eGeom), 'x'+m);
3910 void double_check(t_inputrec *ir, matrix box, gmx_bool bConstr, warninp_t wi)
3914 char warn_buf[STRLEN];
3917 ptr = check_box(ir->ePBC, box);
3920 warning_error(wi, ptr);
3923 if (bConstr && ir->eConstrAlg == econtSHAKE)
3925 if (ir->shake_tol <= 0.0)
3927 sprintf(warn_buf, "ERROR: shake-tol must be > 0 instead of %g\n",
3929 warning_error(wi, warn_buf);
3932 if (IR_TWINRANGE(*ir) && ir->nstlist > 1)
3934 sprintf(warn_buf, "With twin-range cut-off's and SHAKE the virial and the pressure are incorrect.");
3935 if (ir->epc == epcNO)
3937 warning(wi, warn_buf);
3941 warning_error(wi, warn_buf);
3946 if ( (ir->eConstrAlg == econtLINCS) && bConstr)
3948 /* If we have Lincs constraints: */
3949 if (ir->eI == eiMD && ir->etc == etcNO &&
3950 ir->eConstrAlg == econtLINCS && ir->nLincsIter == 1)
3952 sprintf(warn_buf, "For energy conservation with LINCS, lincs_iter should be 2 or larger.\n");
3953 warning_note(wi, warn_buf);
3956 if ((ir->eI == eiCG || ir->eI == eiLBFGS) && (ir->nProjOrder < 8))
3958 sprintf(warn_buf, "For accurate %s with LINCS constraints, lincs-order should be 8 or more.", ei_names[ir->eI]);
3959 warning_note(wi, warn_buf);
3961 if (ir->epc == epcMTTK)
3963 warning_error(wi, "MTTK not compatible with lincs -- use shake instead.");
3967 if (ir->LincsWarnAngle > 90.0)
3969 sprintf(warn_buf, "lincs-warnangle can not be larger than 90 degrees, setting it to 90.\n");
3970 warning(wi, warn_buf);
3971 ir->LincsWarnAngle = 90.0;
3974 if (ir->ePBC != epbcNONE)
3976 if (ir->nstlist == 0)
3978 warning(wi, "With nstlist=0 atoms are only put into the box at step 0, therefore drifting atoms might cause the simulation to crash.");
3980 bTWIN = (ir->rlistlong > ir->rlist);
3981 if (ir->ns_type == ensGRID)
3983 if (sqr(ir->rlistlong) >= max_cutoff2(ir->ePBC, box))
3985 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",
3986 bTWIN ? (ir->rcoulomb == ir->rlistlong ? "rcoulomb" : "rvdw") : "rlist");
3987 warning_error(wi, warn_buf);
3992 min_size = min(box[XX][XX], min(box[YY][YY], box[ZZ][ZZ]));
3993 if (2*ir->rlistlong >= min_size)
3995 sprintf(warn_buf, "ERROR: One of the box lengths is smaller than twice the cut-off length. Increase the box size or decrease rlist.");
3996 warning_error(wi, warn_buf);
3999 fprintf(stderr, "Grid search might allow larger cut-off's than simple search with triclinic boxes.");
4006 void check_chargegroup_radii(const gmx_mtop_t *mtop, const t_inputrec *ir,
4010 real rvdw1, rvdw2, rcoul1, rcoul2;
4011 char warn_buf[STRLEN];
4013 calc_chargegroup_radii(mtop, x, &rvdw1, &rvdw2, &rcoul1, &rcoul2);
4017 printf("Largest charge group radii for Van der Waals: %5.3f, %5.3f nm\n",
4022 printf("Largest charge group radii for Coulomb: %5.3f, %5.3f nm\n",
4028 if (rvdw1 + rvdw2 > ir->rlist ||
4029 rcoul1 + rcoul2 > ir->rlist)
4032 "The sum of the two largest charge group radii (%f) "
4033 "is larger than rlist (%f)\n",
4034 max(rvdw1+rvdw2, rcoul1+rcoul2), ir->rlist);
4035 warning(wi, warn_buf);
4039 /* Here we do not use the zero at cut-off macro,
4040 * since user defined interactions might purposely
4041 * not be zero at the cut-off.
4043 if ((EVDW_IS_ZERO_AT_CUTOFF(ir->vdwtype) ||
4044 ir->vdw_modifier != eintmodNONE) &&
4045 rvdw1 + rvdw2 > ir->rlistlong - ir->rvdw)
4047 sprintf(warn_buf, "The sum of the two largest charge group "
4048 "radii (%f) is larger than %s (%f) - rvdw (%f).\n"
4049 "With exact cut-offs, better performance can be "
4050 "obtained with cutoff-scheme = %s, because it "
4051 "does not use charge groups at all.",
4053 ir->rlistlong > ir->rlist ? "rlistlong" : "rlist",
4054 ir->rlistlong, ir->rvdw,
4055 ecutscheme_names[ecutsVERLET]);
4058 warning(wi, warn_buf);
4062 warning_note(wi, warn_buf);
4065 if ((EEL_IS_ZERO_AT_CUTOFF(ir->coulombtype) ||
4066 ir->coulomb_modifier != eintmodNONE) &&
4067 rcoul1 + rcoul2 > ir->rlistlong - ir->rcoulomb)
4069 sprintf(warn_buf, "The sum of the two largest charge group radii (%f) is larger than %s (%f) - rcoulomb (%f).\n"
4070 "With exact cut-offs, better performance can be obtained with cutoff-scheme = %s, because it does not use charge groups at all.",
4072 ir->rlistlong > ir->rlist ? "rlistlong" : "rlist",
4073 ir->rlistlong, ir->rcoulomb,
4074 ecutscheme_names[ecutsVERLET]);
4077 warning(wi, warn_buf);
4081 warning_note(wi, warn_buf);