2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
49 #include "gmx_fatal.h"
62 #include "mtop_util.h"
63 #include "chargegroup.h"
69 /* Resource parameters
70 * Do not change any of these until you read the instruction
71 * in readinp.h. Some cpp's do not take spaces after the backslash
72 * (like the c-shell), which will give you a very weird compiler
76 typedef struct t_inputrec_strings
78 char tcgrps[STRLEN], tau_t[STRLEN], ref_t[STRLEN],
79 acc[STRLEN], accgrps[STRLEN], freeze[STRLEN], frdim[STRLEN],
80 energy[STRLEN], user1[STRLEN], user2[STRLEN], vcm[STRLEN], xtc_grps[STRLEN],
81 couple_moltype[STRLEN], orirefitgrp[STRLEN], egptable[STRLEN], egpexcl[STRLEN],
82 wall_atomtype[STRLEN], wall_density[STRLEN], deform[STRLEN], QMMM[STRLEN];
83 char fep_lambda[efptNR][STRLEN];
84 char lambda_weights[STRLEN];
87 char anneal[STRLEN], anneal_npoints[STRLEN],
88 anneal_time[STRLEN], anneal_temp[STRLEN];
89 char QMmethod[STRLEN], QMbasis[STRLEN], QMcharge[STRLEN], QMmult[STRLEN],
90 bSH[STRLEN], CASorbitals[STRLEN], CASelectrons[STRLEN], SAon[STRLEN],
91 SAoff[STRLEN], SAsteps[STRLEN], bTS[STRLEN], bOPT[STRLEN];
92 char efield_x[STRLEN], efield_xt[STRLEN], efield_y[STRLEN],
93 efield_yt[STRLEN], efield_z[STRLEN], efield_zt[STRLEN];
95 } gmx_inputrec_strings;
97 static gmx_inputrec_strings *is = NULL;
99 void init_inputrec_strings()
103 gmx_incons("Attempted to call init_inputrec_strings before calling done_inputrec_strings. Only one inputrec (i.e. .mdp file) can be parsed at a time.");
108 void done_inputrec_strings()
115 egrptpALL, /* All particles have to be a member of a group. */
116 egrptpALL_GENREST, /* A rest group with name is generated for particles *
117 * that are not part of any group. */
118 egrptpPART, /* As egrptpALL_GENREST, but no name is generated *
119 * for the rest group. */
120 egrptpONE /* Merge all selected groups into one group, *
121 * make a rest group for the remaining particles. */
124 static const char *constraints[eshNR+1] = {
125 "none", "h-bonds", "all-bonds", "h-angles", "all-angles", NULL
128 static const char *couple_lam[ecouplamNR+1] = {
129 "vdw-q", "vdw", "q", "none", NULL
132 void init_ir(t_inputrec *ir, t_gromppopts *opts)
134 snew(opts->include, STRLEN);
135 snew(opts->define, STRLEN);
136 snew(ir->fepvals, 1);
137 snew(ir->expandedvals, 1);
138 snew(ir->simtempvals, 1);
141 static void GetSimTemps(int ntemps, t_simtemp *simtemp, double *temperature_lambdas)
146 for (i = 0; i < ntemps; i++)
148 /* simple linear scaling -- allows more control */
149 if (simtemp->eSimTempScale == esimtempLINEAR)
151 simtemp->temperatures[i] = simtemp->simtemp_low + (simtemp->simtemp_high-simtemp->simtemp_low)*temperature_lambdas[i];
153 else if (simtemp->eSimTempScale == esimtempGEOMETRIC) /* should give roughly equal acceptance for constant heat capacity . . . */
155 simtemp->temperatures[i] = simtemp->simtemp_low * pow(simtemp->simtemp_high/simtemp->simtemp_low, (1.0*i)/(ntemps-1));
157 else if (simtemp->eSimTempScale == esimtempEXPONENTIAL)
159 simtemp->temperatures[i] = simtemp->simtemp_low + (simtemp->simtemp_high-simtemp->simtemp_low)*((exp(temperature_lambdas[i])-1)/(exp(1.0)-1));
164 sprintf(errorstr, "eSimTempScale=%d not defined", simtemp->eSimTempScale);
165 gmx_fatal(FARGS, errorstr);
172 static void _low_check(gmx_bool b, char *s, warninp_t wi)
176 warning_error(wi, s);
180 static void check_nst(const char *desc_nst, int nst,
181 const char *desc_p, int *p,
186 if (*p > 0 && *p % nst != 0)
188 /* Round up to the next multiple of nst */
189 *p = ((*p)/nst + 1)*nst;
190 sprintf(buf, "%s should be a multiple of %s, changing %s to %d\n",
191 desc_p, desc_nst, desc_p, *p);
196 static gmx_bool ir_NVE(const t_inputrec *ir)
198 return ((ir->eI == eiMD || EI_VV(ir->eI)) && ir->etc == etcNO);
201 static int lcd(int n1, int n2)
206 for (i = 2; (i <= n1 && i <= n2); i++)
208 if (n1 % i == 0 && n2 % i == 0)
217 static void process_interaction_modifier(const t_inputrec *ir, int *eintmod)
219 if (*eintmod == eintmodPOTSHIFT_VERLET)
221 if (ir->cutoff_scheme == ecutsVERLET)
223 *eintmod = eintmodPOTSHIFT;
227 *eintmod = eintmodNONE;
232 void check_ir(const char *mdparin, t_inputrec *ir, t_gromppopts *opts,
234 /* Check internal consistency */
236 /* Strange macro: first one fills the err_buf, and then one can check
237 * the condition, which will print the message and increase the error
240 #define CHECK(b) _low_check(b, err_buf, wi)
241 char err_buf[256], warn_buf[STRLEN];
247 t_lambda *fep = ir->fepvals;
248 t_expanded *expand = ir->expandedvals;
250 set_warning_line(wi, mdparin, -1);
252 /* BASIC CUT-OFF STUFF */
253 if (ir->rcoulomb < 0)
255 warning_error(wi, "rcoulomb should be >= 0");
259 warning_error(wi, "rvdw should be >= 0");
262 !(ir->cutoff_scheme == ecutsVERLET && ir->verletbuf_tol > 0))
264 warning_error(wi, "rlist should be >= 0");
267 process_interaction_modifier(ir, &ir->coulomb_modifier);
268 process_interaction_modifier(ir, &ir->vdw_modifier);
270 if (ir->cutoff_scheme == ecutsGROUP)
272 /* BASIC CUT-OFF STUFF */
273 if (ir->rlist == 0 ||
274 !((EEL_MIGHT_BE_ZERO_AT_CUTOFF(ir->coulombtype) && ir->rcoulomb > ir->rlist) ||
275 (EVDW_MIGHT_BE_ZERO_AT_CUTOFF(ir->vdwtype) && ir->rvdw > ir->rlist)))
277 /* No switched potential and/or no twin-range:
278 * we can set the long-range cut-off to the maximum of the other cut-offs.
280 ir->rlistlong = max_cutoff(ir->rlist, max_cutoff(ir->rvdw, ir->rcoulomb));
282 else if (ir->rlistlong < 0)
284 ir->rlistlong = max_cutoff(ir->rlist, max_cutoff(ir->rvdw, ir->rcoulomb));
285 sprintf(warn_buf, "rlistlong was not set, setting it to %g (no buffer)",
287 warning(wi, warn_buf);
289 if (ir->rlistlong == 0 && ir->ePBC != epbcNONE)
291 warning_error(wi, "Can not have an infinite cut-off with PBC");
293 if (ir->rlistlong > 0 && (ir->rlist == 0 || ir->rlistlong < ir->rlist))
295 warning_error(wi, "rlistlong can not be shorter than rlist");
297 if (IR_TWINRANGE(*ir) && ir->nstlist <= 0)
299 warning_error(wi, "Can not have nstlist<=0 with twin-range interactions");
303 if (ir->rlistlong == ir->rlist)
307 else if (ir->rlistlong > ir->rlist && ir->nstcalclr == 0)
309 warning_error(wi, "With different cutoffs for electrostatics and VdW, nstcalclr must be -1 or a positive number");
312 if (ir->cutoff_scheme == ecutsVERLET)
316 /* Normal Verlet type neighbor-list, currently only limited feature support */
317 if (inputrec2nboundeddim(ir) < 3)
319 warning_error(wi, "With Verlet lists only full pbc or pbc=xy with walls is supported");
321 if (ir->rcoulomb != ir->rvdw)
323 warning_error(wi, "With Verlet lists rcoulomb!=rvdw is not supported");
325 if (ir->vdwtype != evdwCUT)
327 warning_error(wi, "With Verlet lists only cut-off LJ interactions are supported");
329 if (!(ir->coulombtype == eelCUT ||
330 (EEL_RF(ir->coulombtype) && ir->coulombtype != eelRF_NEC) ||
331 EEL_PME(ir->coulombtype) || ir->coulombtype == eelEWALD))
333 warning_error(wi, "With Verlet lists only cut-off, reaction-field, PME and Ewald electrostatics are supported");
336 if (ir->nstlist <= 0)
338 warning_error(wi, "With Verlet lists nstlist should be larger than 0");
341 if (ir->nstlist < 10)
343 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.");
346 rc_max = max(ir->rvdw, ir->rcoulomb);
348 if (ir->verletbuf_tol <= 0)
350 if (ir->verletbuf_tol == 0)
352 warning_error(wi, "Can not have Verlet buffer tolerance of exactly 0");
355 if (ir->rlist < rc_max)
357 warning_error(wi, "With verlet lists rlist can not be smaller than rvdw or rcoulomb");
360 if (ir->rlist == rc_max && ir->nstlist > 1)
362 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.");
367 if (ir->rlist > rc_max)
369 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.");
372 if (ir->nstlist == 1)
374 /* No buffer required */
379 if (EI_DYNAMICS(ir->eI))
381 if (EI_MD(ir->eI) && ir->etc == etcNO)
383 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.");
386 if (inputrec2nboundeddim(ir) < 3)
388 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.");
390 /* Set rlist temporarily so we can continue processing */
395 /* Set the buffer to 5% of the cut-off */
396 ir->rlist = 1.05*rc_max;
401 /* No twin-range calculations with Verlet lists */
402 ir->rlistlong = ir->rlist;
405 if (ir->nstcalclr == -1)
407 /* if rlist=rlistlong, this will later be changed to nstcalclr=0 */
408 ir->nstcalclr = ir->nstlist;
410 else if (ir->nstcalclr > 0)
412 if (ir->nstlist > 0 && (ir->nstlist % ir->nstcalclr != 0))
414 warning_error(wi, "nstlist must be evenly divisible by nstcalclr. Use nstcalclr = -1 to automatically follow nstlist");
417 else if (ir->nstcalclr < -1)
419 warning_error(wi, "nstcalclr must be a positive number (divisor of nstcalclr), or -1 to follow nstlist.");
422 if (EEL_PME(ir->coulombtype) && ir->rcoulomb > ir->rvdw && ir->nstcalclr > 1)
424 warning_error(wi, "When used with PME, the long-range component of twin-range interactions must be updated every step (nstcalclr)");
427 /* GENERAL INTEGRATOR STUFF */
428 if (!(ir->eI == eiMD || EI_VV(ir->eI)))
432 if (ir->eI == eiVVAK)
434 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]);
435 warning_note(wi, warn_buf);
437 if (!EI_DYNAMICS(ir->eI))
441 if (EI_DYNAMICS(ir->eI))
443 if (ir->nstcalcenergy < 0)
445 ir->nstcalcenergy = ir_optimal_nstcalcenergy(ir);
446 if (ir->nstenergy != 0 && ir->nstenergy < ir->nstcalcenergy)
448 /* nstcalcenergy larger than nstener does not make sense.
449 * We ideally want nstcalcenergy=nstener.
453 ir->nstcalcenergy = lcd(ir->nstenergy, ir->nstlist);
457 ir->nstcalcenergy = ir->nstenergy;
461 else if ( (ir->nstenergy > 0 && ir->nstcalcenergy > ir->nstenergy) ||
462 (ir->efep != efepNO && ir->fepvals->nstdhdl > 0 &&
463 (ir->nstcalcenergy > ir->fepvals->nstdhdl) ) )
466 const char *nsten = "nstenergy";
467 const char *nstdh = "nstdhdl";
468 const char *min_name = nsten;
469 int min_nst = ir->nstenergy;
471 /* find the smallest of ( nstenergy, nstdhdl ) */
472 if (ir->efep != efepNO && ir->fepvals->nstdhdl > 0 &&
473 (ir->fepvals->nstdhdl < ir->nstenergy) )
475 min_nst = ir->fepvals->nstdhdl;
478 /* If the user sets nstenergy small, we should respect that */
480 "Setting nstcalcenergy (%d) equal to %s (%d)",
481 ir->nstcalcenergy, min_name, min_nst);
482 warning_note(wi, warn_buf);
483 ir->nstcalcenergy = min_nst;
486 if (ir->epc != epcNO)
488 if (ir->nstpcouple < 0)
490 ir->nstpcouple = ir_optimal_nstpcouple(ir);
493 if (IR_TWINRANGE(*ir))
495 check_nst("nstlist", ir->nstlist,
496 "nstcalcenergy", &ir->nstcalcenergy, wi);
497 if (ir->epc != epcNO)
499 check_nst("nstlist", ir->nstlist,
500 "nstpcouple", &ir->nstpcouple, wi);
504 if (ir->nstcalcenergy > 0)
506 if (ir->efep != efepNO)
508 /* nstdhdl should be a multiple of nstcalcenergy */
509 check_nst("nstcalcenergy", ir->nstcalcenergy,
510 "nstdhdl", &ir->fepvals->nstdhdl, wi);
511 /* nstexpanded should be a multiple of nstcalcenergy */
512 check_nst("nstcalcenergy", ir->nstcalcenergy,
513 "nstexpanded", &ir->expandedvals->nstexpanded, wi);
515 /* for storing exact averages nstenergy should be
516 * a multiple of nstcalcenergy
518 check_nst("nstcalcenergy", ir->nstcalcenergy,
519 "nstenergy", &ir->nstenergy, wi);
524 if ((EI_SD(ir->eI) || ir->eI == eiBD) &&
525 ir->bContinuation && ir->ld_seed != -1)
527 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)");
533 sprintf(err_buf, "TPI only works with pbc = %s", epbc_names[epbcXYZ]);
534 CHECK(ir->ePBC != epbcXYZ);
535 sprintf(err_buf, "TPI only works with ns = %s", ens_names[ensGRID]);
536 CHECK(ir->ns_type != ensGRID);
537 sprintf(err_buf, "with TPI nstlist should be larger than zero");
538 CHECK(ir->nstlist <= 0);
539 sprintf(err_buf, "TPI does not work with full electrostatics other than PME");
540 CHECK(EEL_FULL(ir->coulombtype) && !EEL_PME(ir->coulombtype));
544 if ( (opts->nshake > 0) && (opts->bMorse) )
547 "Using morse bond-potentials while constraining bonds is useless");
548 warning(wi, warn_buf);
551 if ((EI_SD(ir->eI) || ir->eI == eiBD) &&
552 ir->bContinuation && ir->ld_seed != -1)
554 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)");
556 /* verify simulated tempering options */
560 gmx_bool bAllTempZero = TRUE;
561 for (i = 0; i < fep->n_lambda; i++)
563 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]);
564 CHECK((fep->all_lambda[efptTEMPERATURE][i] < 0) || (fep->all_lambda[efptTEMPERATURE][i] > 1));
565 if (fep->all_lambda[efptTEMPERATURE][i] > 0)
567 bAllTempZero = FALSE;
570 sprintf(err_buf, "if simulated tempering is on, temperature-lambdas may not be all zero");
571 CHECK(bAllTempZero == TRUE);
573 sprintf(err_buf, "Simulated tempering is currently only compatible with md-vv");
574 CHECK(ir->eI != eiVV);
576 /* check compatability of the temperature coupling with simulated tempering */
578 if (ir->etc == etcNOSEHOOVER)
580 sprintf(warn_buf, "Nose-Hoover based temperature control such as [%s] my not be entirelyconsistent with simulated tempering", etcoupl_names[ir->etc]);
581 warning_note(wi, warn_buf);
584 /* check that the temperatures make sense */
586 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);
587 CHECK(ir->simtempvals->simtemp_high <= ir->simtempvals->simtemp_low);
589 sprintf(err_buf, "Higher simulated tempering temperature (%g) must be >= zero", ir->simtempvals->simtemp_high);
590 CHECK(ir->simtempvals->simtemp_high <= 0);
592 sprintf(err_buf, "Lower simulated tempering temperature (%g) must be >= zero", ir->simtempvals->simtemp_low);
593 CHECK(ir->simtempvals->simtemp_low <= 0);
596 /* verify free energy options */
598 if (ir->efep != efepNO)
601 sprintf(err_buf, "The soft-core power is %d and can only be 1 or 2",
603 CHECK(fep->sc_alpha != 0 && fep->sc_power != 1 && fep->sc_power != 2);
605 sprintf(err_buf, "The soft-core sc-r-power is %d and can only be 6 or 48",
606 (int)fep->sc_r_power);
607 CHECK(fep->sc_alpha != 0 && fep->sc_r_power != 6.0 && fep->sc_r_power != 48.0);
609 sprintf(err_buf, "Can't use postive delta-lambda (%g) if initial state/lambda does not start at zero", fep->delta_lambda);
610 CHECK(fep->delta_lambda > 0 && ((fep->init_fep_state > 0) || (fep->init_lambda > 0)));
612 sprintf(err_buf, "Can't use postive delta-lambda (%g) with expanded ensemble simulations", fep->delta_lambda);
613 CHECK(fep->delta_lambda > 0 && (ir->efep == efepEXPANDED));
615 sprintf(err_buf, "Can only use expanded ensemble with md-vv for now; should be supported for other integrators in 5.0");
616 CHECK(!(EI_VV(ir->eI)) && (ir->efep == efepEXPANDED));
618 sprintf(err_buf, "Free-energy not implemented for Ewald");
619 CHECK(ir->coulombtype == eelEWALD);
621 /* check validty of lambda inputs */
622 if (fep->n_lambda == 0)
624 /* Clear output in case of no states:*/
625 sprintf(err_buf, "init-lambda-state set to %d: no lambda states are defined.", fep->init_fep_state);
626 CHECK((fep->init_fep_state >= 0) && (fep->n_lambda == 0));
630 sprintf(err_buf, "initial thermodynamic state %d does not exist, only goes to %d", fep->init_fep_state, fep->n_lambda-1);
631 CHECK((fep->init_fep_state >= fep->n_lambda));
634 sprintf(err_buf, "Lambda state must be set, either with init-lambda-state or with init-lambda");
635 CHECK((fep->init_fep_state < 0) && (fep->init_lambda < 0));
637 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",
638 fep->init_lambda, fep->init_fep_state);
639 CHECK((fep->init_fep_state >= 0) && (fep->init_lambda >= 0));
643 if ((fep->init_lambda >= 0) && (fep->delta_lambda == 0))
647 for (i = 0; i < efptNR; i++)
649 if (fep->separate_dvdl[i])
654 if (n_lambda_terms > 1)
656 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.");
657 warning(wi, warn_buf);
660 if (n_lambda_terms < 2 && fep->n_lambda > 0)
663 "init-lambda is deprecated for setting lambda state (except for slow growth). Use init-lambda-state instead.");
667 for (j = 0; j < efptNR; j++)
669 for (i = 0; i < fep->n_lambda; i++)
671 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]);
672 CHECK((fep->all_lambda[j][i] < 0) || (fep->all_lambda[j][i] > 1));
676 if ((fep->sc_alpha > 0) && (!fep->bScCoul))
678 for (i = 0; i < fep->n_lambda; i++)
680 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],
681 fep->all_lambda[efptCOUL][i]);
682 CHECK((fep->sc_alpha > 0) &&
683 (((fep->all_lambda[efptCOUL][i] > 0.0) &&
684 (fep->all_lambda[efptCOUL][i] < 1.0)) &&
685 ((fep->all_lambda[efptVDW][i] > 0.0) &&
686 (fep->all_lambda[efptVDW][i] < 1.0))));
690 if ((fep->bScCoul) && (EEL_PME(ir->coulombtype)))
692 real sigma, lambda, r_sc;
695 /* Maximum estimate for A and B charges equal with lambda power 1 */
697 r_sc = pow(lambda*fep->sc_alpha*pow(sigma/ir->rcoulomb, fep->sc_r_power) + 1.0, 1.0/fep->sc_r_power);
698 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.",
700 sigma, lambda, r_sc - 1.0, ir->ewald_rtol);
701 warning_note(wi, warn_buf);
704 /* Free Energy Checks -- In an ideal world, slow growth and FEP would
705 be treated differently, but that's the next step */
707 for (i = 0; i < efptNR; i++)
709 for (j = 0; j < fep->n_lambda; j++)
711 sprintf(err_buf, "%s[%d] must be between 0 and 1", efpt_names[i], j);
712 CHECK((fep->all_lambda[i][j] < 0) || (fep->all_lambda[i][j] > 1));
717 if ((ir->bSimTemp) || (ir->efep == efepEXPANDED))
720 expand = ir->expandedvals;
722 /* checking equilibration of weights inputs for validity */
724 sprintf(err_buf, "weight-equil-number-all-lambda (%d) is ignored if lmc-weights-equil is not equal to %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) is ignored if lmc-weights-equil is not equal to %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) is ignored if lmc-weights-equil is not equal to %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 (%d) is ignored if lmc-weights-equil is not equal to %s",
737 expand->equil_samples, elmceq_names[elmceqWLDELTA]);
738 CHECK((expand->equil_wl_delta > 0) && (expand->elmceq != elmceqWLDELTA));
740 sprintf(err_buf, "weight-equil-count-ratio (%f) is ignored if lmc-weights-equil is not equal to %s",
741 expand->equil_ratio, elmceq_names[elmceqRATIO]);
742 CHECK((expand->equil_ratio > 0) && (expand->elmceq != elmceqRATIO));
744 sprintf(err_buf, "weight-equil-number-all-lambda (%d) must be a positive integer if lmc-weights-equil=%s",
745 expand->equil_n_at_lam, elmceq_names[elmceqNUMATLAM]);
746 CHECK((expand->equil_n_at_lam <= 0) && (expand->elmceq == elmceqNUMATLAM));
748 sprintf(err_buf, "weight-equil-number-samples (%d) must be a positive integer if lmc-weights-equil=%s",
749 expand->equil_samples, elmceq_names[elmceqSAMPLES]);
750 CHECK((expand->equil_samples <= 0) && (expand->elmceq == elmceqSAMPLES));
752 sprintf(err_buf, "weight-equil-number-steps (%d) must be a positive integer if lmc-weights-equil=%s",
753 expand->equil_steps, elmceq_names[elmceqSTEPS]);
754 CHECK((expand->equil_steps <= 0) && (expand->elmceq == elmceqSTEPS));
756 sprintf(err_buf, "weight-equil-wl-delta (%f) must be > 0 if lmc-weights-equil=%s",
757 expand->equil_wl_delta, elmceq_names[elmceqWLDELTA]);
758 CHECK((expand->equil_wl_delta <= 0) && (expand->elmceq == elmceqWLDELTA));
760 sprintf(err_buf, "weight-equil-count-ratio (%f) must be > 0 if lmc-weights-equil=%s",
761 expand->equil_ratio, elmceq_names[elmceqRATIO]);
762 CHECK((expand->equil_ratio <= 0) && (expand->elmceq == elmceqRATIO));
764 sprintf(err_buf, "lmc-weights-equil=%s only possible when lmc-stats = %s or lmc-stats %s",
765 elmceq_names[elmceqWLDELTA], elamstats_names[elamstatsWL], elamstats_names[elamstatsWWL]);
766 CHECK((expand->elmceq == elmceqWLDELTA) && (!EWL(expand->elamstats)));
768 sprintf(err_buf, "lmc-repeats (%d) must be greater than 0", expand->lmc_repeats);
769 CHECK((expand->lmc_repeats <= 0));
770 sprintf(err_buf, "minimum-var-min (%d) must be greater than 0", expand->minvarmin);
771 CHECK((expand->minvarmin <= 0));
772 sprintf(err_buf, "weight-c-range (%d) must be greater or equal to 0", expand->c_range);
773 CHECK((expand->c_range < 0));
774 sprintf(err_buf, "init-lambda-state (%d) must be zero if lmc-forced-nstart (%d)> 0 and lmc-move != 'no'",
775 fep->init_fep_state, expand->lmc_forced_nstart);
776 CHECK((fep->init_fep_state != 0) && (expand->lmc_forced_nstart > 0) && (expand->elmcmove != elmcmoveNO));
777 sprintf(err_buf, "lmc-forced-nstart (%d) must not be negative", expand->lmc_forced_nstart);
778 CHECK((expand->lmc_forced_nstart < 0));
779 sprintf(err_buf, "init-lambda-state (%d) must be in the interval [0,number of lambdas)", fep->init_fep_state);
780 CHECK((fep->init_fep_state < 0) || (fep->init_fep_state >= fep->n_lambda));
782 sprintf(err_buf, "init-wl-delta (%f) must be greater than or equal to 0", expand->init_wl_delta);
783 CHECK((expand->init_wl_delta < 0));
784 sprintf(err_buf, "wl-ratio (%f) must be between 0 and 1", expand->wl_ratio);
785 CHECK((expand->wl_ratio <= 0) || (expand->wl_ratio >= 1));
786 sprintf(err_buf, "wl-scale (%f) must be between 0 and 1", expand->wl_scale);
787 CHECK((expand->wl_scale <= 0) || (expand->wl_scale >= 1));
789 /* if there is no temperature control, we need to specify an MC temperature */
790 sprintf(err_buf, "If there is no temperature control, and lmc-mcmove!= 'no',mc_temperature must be set to a positive number");
791 if (expand->nstTij > 0)
793 sprintf(err_buf, "nst-transition-matrix (%d) must be an integer multiple of nstlog (%d)",
794 expand->nstTij, ir->nstlog);
795 CHECK((mod(expand->nstTij, ir->nstlog) != 0));
800 sprintf(err_buf, "walls only work with pbc=%s", epbc_names[epbcXY]);
801 CHECK(ir->nwall && ir->ePBC != epbcXY);
804 if (ir->ePBC != epbcXYZ && ir->nwall != 2)
806 if (ir->ePBC == epbcNONE)
808 if (ir->epc != epcNO)
810 warning(wi, "Turning off pressure coupling for vacuum system");
816 sprintf(err_buf, "Can not have pressure coupling with pbc=%s",
817 epbc_names[ir->ePBC]);
818 CHECK(ir->epc != epcNO);
820 sprintf(err_buf, "Can not have Ewald with pbc=%s", epbc_names[ir->ePBC]);
821 CHECK(EEL_FULL(ir->coulombtype));
823 sprintf(err_buf, "Can not have dispersion correction with pbc=%s",
824 epbc_names[ir->ePBC]);
825 CHECK(ir->eDispCorr != edispcNO);
828 if (ir->rlist == 0.0)
830 sprintf(err_buf, "can only have neighborlist cut-off zero (=infinite)\n"
831 "with coulombtype = %s or coulombtype = %s\n"
832 "without periodic boundary conditions (pbc = %s) and\n"
833 "rcoulomb and rvdw set to zero",
834 eel_names[eelCUT], eel_names[eelUSER], epbc_names[epbcNONE]);
835 CHECK(((ir->coulombtype != eelCUT) && (ir->coulombtype != eelUSER)) ||
836 (ir->ePBC != epbcNONE) ||
837 (ir->rcoulomb != 0.0) || (ir->rvdw != 0.0));
841 warning_error(wi, "Can not have heuristic neighborlist updates without cut-off");
845 warning_note(wi, "Simulating without cut-offs is usually (slightly) faster with nstlist=0, nstype=simple and particle decomposition");
850 if (ir->nstcomm == 0)
852 ir->comm_mode = ecmNO;
854 if (ir->comm_mode != ecmNO)
858 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");
859 ir->nstcomm = abs(ir->nstcomm);
862 if (ir->nstcalcenergy > 0 && ir->nstcomm < ir->nstcalcenergy)
864 warning_note(wi, "nstcomm < nstcalcenergy defeats the purpose of nstcalcenergy, setting nstcomm to nstcalcenergy");
865 ir->nstcomm = ir->nstcalcenergy;
868 if (ir->comm_mode == ecmANGULAR)
870 sprintf(err_buf, "Can not remove the rotation around the center of mass with periodic molecules");
871 CHECK(ir->bPeriodicMols);
872 if (ir->ePBC != epbcNONE)
874 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).");
879 if (EI_STATE_VELOCITY(ir->eI) && ir->ePBC == epbcNONE && ir->comm_mode != ecmANGULAR)
881 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.");
884 sprintf(err_buf, "Twin-range neighbour searching (NS) with simple NS"
885 " algorithm not implemented");
886 CHECK(((ir->rcoulomb > ir->rlist) || (ir->rvdw > ir->rlist))
887 && (ir->ns_type == ensSIMPLE));
889 /* TEMPERATURE COUPLING */
890 if (ir->etc == etcYES)
892 ir->etc = etcBERENDSEN;
893 warning_note(wi, "Old option for temperature coupling given: "
894 "changing \"yes\" to \"Berendsen\"\n");
897 if ((ir->etc == etcNOSEHOOVER) || (ir->epc == epcMTTK))
899 if (ir->opts.nhchainlength < 1)
901 sprintf(warn_buf, "number of Nose-Hoover chains (currently %d) cannot be less than 1,reset to 1\n", ir->opts.nhchainlength);
902 ir->opts.nhchainlength = 1;
903 warning(wi, warn_buf);
906 if (ir->etc == etcNOSEHOOVER && !EI_VV(ir->eI) && ir->opts.nhchainlength > 1)
908 warning_note(wi, "leapfrog does not yet support Nose-Hoover chains, nhchainlength reset to 1");
909 ir->opts.nhchainlength = 1;
914 ir->opts.nhchainlength = 0;
917 if (ir->eI == eiVVAK)
919 sprintf(err_buf, "%s implemented primarily for validation, and requires nsttcouple = 1 and nstpcouple = 1.",
921 CHECK((ir->nsttcouple != 1) || (ir->nstpcouple != 1));
924 if (ETC_ANDERSEN(ir->etc))
926 sprintf(err_buf, "%s temperature control not supported for integrator %s.", etcoupl_names[ir->etc], ei_names[ir->eI]);
927 CHECK(!(EI_VV(ir->eI)));
929 for (i = 0; i < ir->opts.ngtc; i++)
931 sprintf(err_buf, "all tau_t must currently be equal using Andersen temperature control, violated for group %d", i);
932 CHECK(ir->opts.tau_t[0] != ir->opts.tau_t[i]);
933 sprintf(err_buf, "all tau_t must be postive using Andersen temperature control, tau_t[%d]=%10.6f",
934 i, ir->opts.tau_t[i]);
935 CHECK(ir->opts.tau_t[i] < 0);
937 if (ir->nstcomm > 0 && (ir->etc == etcANDERSEN))
939 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]);
940 warning_note(wi, warn_buf);
943 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]);
944 CHECK(ir->nstcomm > 1 && (ir->etc == etcANDERSEN));
946 for (i = 0; i < ir->opts.ngtc; i++)
948 int nsteps = (int)(ir->opts.tau_t[i]/ir->delta_t);
949 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);
950 CHECK((nsteps % ir->nstcomm) && (ir->etc == etcANDERSENMASSIVE));
953 if (ir->etc == etcBERENDSEN)
955 sprintf(warn_buf, "The %s thermostat does not generate the correct kinetic energy distribution. You might want to consider using the %s thermostat.",
956 ETCOUPLTYPE(ir->etc), ETCOUPLTYPE(etcVRESCALE));
957 warning_note(wi, warn_buf);
960 if ((ir->etc == etcNOSEHOOVER || ETC_ANDERSEN(ir->etc))
961 && ir->epc == epcBERENDSEN)
963 sprintf(warn_buf, "Using Berendsen pressure coupling invalidates the "
964 "true ensemble for the thermostat");
965 warning(wi, warn_buf);
968 /* PRESSURE COUPLING */
969 if (ir->epc == epcISOTROPIC)
971 ir->epc = epcBERENDSEN;
972 warning_note(wi, "Old option for pressure coupling given: "
973 "changing \"Isotropic\" to \"Berendsen\"\n");
976 if (ir->epc != epcNO)
978 dt_pcoupl = ir->nstpcouple*ir->delta_t;
980 sprintf(err_buf, "tau-p must be > 0 instead of %g\n", ir->tau_p);
981 CHECK(ir->tau_p <= 0);
983 if (ir->tau_p/dt_pcoupl < pcouple_min_integration_steps(ir->epc))
985 sprintf(warn_buf, "For proper integration of the %s barostat, tau-p (%g) should be at least %d times larger than nstpcouple*dt (%g)",
986 EPCOUPLTYPE(ir->epc), ir->tau_p, pcouple_min_integration_steps(ir->epc), dt_pcoupl);
987 warning(wi, warn_buf);
990 sprintf(err_buf, "compressibility must be > 0 when using pressure"
991 " coupling %s\n", EPCOUPLTYPE(ir->epc));
992 CHECK(ir->compress[XX][XX] < 0 || ir->compress[YY][YY] < 0 ||
993 ir->compress[ZZ][ZZ] < 0 ||
994 (trace(ir->compress) == 0 && ir->compress[YY][XX] <= 0 &&
995 ir->compress[ZZ][XX] <= 0 && ir->compress[ZZ][YY] <= 0));
997 if (epcPARRINELLORAHMAN == ir->epc && opts->bGenVel)
1000 "You are generating velocities so I am assuming you "
1001 "are equilibrating a system. You are using "
1002 "%s pressure coupling, but this can be "
1003 "unstable for equilibration. If your system crashes, try "
1004 "equilibrating first with Berendsen pressure coupling. If "
1005 "you are not equilibrating the system, you can probably "
1006 "ignore this warning.",
1007 epcoupl_names[ir->epc]);
1008 warning(wi, warn_buf);
1014 if (ir->epc > epcNO)
1016 if ((ir->epc != epcBERENDSEN) && (ir->epc != epcMTTK))
1018 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.");
1023 /* ELECTROSTATICS */
1024 /* More checks are in triple check (grompp.c) */
1026 if (ir->coulombtype == eelSWITCH)
1028 sprintf(warn_buf, "coulombtype = %s is only for testing purposes and can lead to serious "
1029 "artifacts, advice: use coulombtype = %s",
1030 eel_names[ir->coulombtype],
1031 eel_names[eelRF_ZERO]);
1032 warning(wi, warn_buf);
1035 if (ir->epsilon_r != 1 && ir->implicit_solvent == eisGBSA)
1037 sprintf(warn_buf, "epsilon-r = %g with GB implicit solvent, will use this value for inner dielectric", ir->epsilon_r);
1038 warning_note(wi, warn_buf);
1041 if (EEL_RF(ir->coulombtype) && ir->epsilon_rf == 1 && ir->epsilon_r != 1)
1043 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);
1044 warning(wi, warn_buf);
1045 ir->epsilon_rf = ir->epsilon_r;
1046 ir->epsilon_r = 1.0;
1049 if (getenv("GALACTIC_DYNAMICS") == NULL)
1051 sprintf(err_buf, "epsilon-r must be >= 0 instead of %g\n", ir->epsilon_r);
1052 CHECK(ir->epsilon_r < 0);
1055 if (EEL_RF(ir->coulombtype))
1057 /* reaction field (at the cut-off) */
1059 if (ir->coulombtype == eelRF_ZERO)
1061 sprintf(warn_buf, "With coulombtype = %s, epsilon-rf must be 0, assuming you meant epsilon_rf=0",
1062 eel_names[ir->coulombtype]);
1063 CHECK(ir->epsilon_rf != 0);
1064 ir->epsilon_rf = 0.0;
1067 sprintf(err_buf, "epsilon-rf must be >= epsilon-r");
1068 CHECK((ir->epsilon_rf < ir->epsilon_r && ir->epsilon_rf != 0) ||
1069 (ir->epsilon_r == 0));
1070 if (ir->epsilon_rf == ir->epsilon_r)
1072 sprintf(warn_buf, "Using epsilon-rf = epsilon-r with %s does not make sense",
1073 eel_names[ir->coulombtype]);
1074 warning(wi, warn_buf);
1077 /* Allow rlist>rcoulomb for tabulated long range stuff. This just
1078 * means the interaction is zero outside rcoulomb, but it helps to
1079 * provide accurate energy conservation.
1081 if (EEL_MIGHT_BE_ZERO_AT_CUTOFF(ir->coulombtype))
1083 if (EEL_SWITCHED(ir->coulombtype))
1086 "With coulombtype = %s rcoulomb_switch must be < rcoulomb. Or, better: Use the potential modifier options!",
1087 eel_names[ir->coulombtype]);
1088 CHECK(ir->rcoulomb_switch >= ir->rcoulomb);
1091 else if (ir->coulombtype == eelCUT || EEL_RF(ir->coulombtype))
1093 if (ir->cutoff_scheme == ecutsGROUP && ir->coulomb_modifier == eintmodNONE)
1095 sprintf(err_buf, "With coulombtype = %s, rcoulomb should be >= rlist unless you use a potential modifier",
1096 eel_names[ir->coulombtype]);
1097 CHECK(ir->rlist > ir->rcoulomb);
1101 if (ir->coulombtype == eelSWITCH || ir->coulombtype == eelSHIFT ||
1102 ir->vdwtype == evdwSWITCH || ir->vdwtype == evdwSHIFT)
1105 "The switch/shift interaction settings are just for compatibility; you will get better "
1106 "performance from applying potential modifiers to your interactions!\n");
1107 warning_note(wi, warn_buf);
1110 if (ir->coulombtype == eelPMESWITCH)
1112 if (ir->rcoulomb_switch/ir->rcoulomb < 0.9499)
1114 sprintf(warn_buf, "The switching range for %s should be 5%% or less, energy conservation will be good anyhow, since ewald_rtol = %g",
1115 eel_names[ir->coulombtype],
1117 warning(wi, warn_buf);
1121 if (EEL_FULL(ir->coulombtype))
1123 if (ir->coulombtype == eelPMESWITCH || ir->coulombtype == eelPMEUSER ||
1124 ir->coulombtype == eelPMEUSERSWITCH)
1126 sprintf(err_buf, "With coulombtype = %s, rcoulomb must be <= rlist",
1127 eel_names[ir->coulombtype]);
1128 CHECK(ir->rcoulomb > ir->rlist);
1130 else if (ir->cutoff_scheme == ecutsGROUP && ir->coulomb_modifier == eintmodNONE)
1132 if (ir->coulombtype == eelPME || ir->coulombtype == eelP3M_AD)
1135 "With coulombtype = %s (without modifier), rcoulomb must be equal to rlist,\n"
1136 "or rlistlong if nstcalclr=1. For optimal energy conservation,consider using\n"
1137 "a potential modifier.", eel_names[ir->coulombtype]);
1138 if (ir->nstcalclr == 1)
1140 CHECK(ir->rcoulomb != ir->rlist && ir->rcoulomb != ir->rlistlong);
1144 CHECK(ir->rcoulomb != ir->rlist);
1150 if (EVDW_PME(ir->vdwtype))
1152 if (EVDW_MIGHT_BE_ZERO_AT_CUTOFF(ir->vdwtype))
1154 sprintf(err_buf, "With vdwtype = %s, rvdw must be <= rlist",
1155 evdw_names[ir->vdwtype]);
1156 CHECK(ir->rvdw > ir->rlist);
1161 "With vdwtype = %s, rvdw must be equal to rlist\n",
1162 evdw_names[ir->vdwtype]);
1163 CHECK(ir->rvdw != ir->rlist);
1167 if (EEL_PME(ir->coulombtype) || EVDW_PME(ir->vdwtype))
1169 if (ir->pme_order < 3)
1171 warning_error(wi, "pme-order can not be smaller than 3");
1175 if (ir->nwall == 2 && EEL_FULL(ir->coulombtype))
1177 if (ir->ewald_geometry == eewg3D)
1179 sprintf(warn_buf, "With pbc=%s you should use ewald-geometry=%s",
1180 epbc_names[ir->ePBC], eewg_names[eewg3DC]);
1181 warning(wi, warn_buf);
1183 /* This check avoids extra pbc coding for exclusion corrections */
1184 sprintf(err_buf, "wall-ewald-zfac should be >= 2");
1185 CHECK(ir->wall_ewald_zfac < 2);
1188 if (EVDW_SWITCHED(ir->vdwtype))
1190 sprintf(err_buf, "With vdwtype = %s rvdw-switch must be < rvdw. Or, better - use a potential modifier.",
1191 evdw_names[ir->vdwtype]);
1192 CHECK(ir->rvdw_switch >= ir->rvdw);
1194 else if (ir->vdwtype == evdwCUT)
1196 if (ir->cutoff_scheme == ecutsGROUP && ir->vdw_modifier == eintmodNONE)
1198 sprintf(err_buf, "With vdwtype = %s, rvdw must be >= rlist unless you use a potential modifier", evdw_names[ir->vdwtype]);
1199 CHECK(ir->rlist > ir->rvdw);
1202 if (ir->cutoff_scheme == ecutsGROUP)
1204 if (((ir->coulomb_modifier != eintmodNONE && ir->rcoulomb == ir->rlist) ||
1205 (ir->vdw_modifier != eintmodNONE && ir->rvdw == ir->rlist)) &&
1208 warning_note(wi, "With exact cut-offs, rlist should be "
1209 "larger than rcoulomb and rvdw, so that there "
1210 "is a buffer region for particle motion "
1211 "between neighborsearch steps");
1214 if (EEL_IS_ZERO_AT_CUTOFF(ir->coulombtype)
1215 && (ir->rlistlong <= ir->rcoulomb))
1217 sprintf(warn_buf, "For energy conservation with switch/shift potentials, %s should be 0.1 to 0.3 nm larger than rcoulomb.",
1218 IR_TWINRANGE(*ir) ? "rlistlong" : "rlist");
1219 warning_note(wi, warn_buf);
1221 if (EVDW_SWITCHED(ir->vdwtype) && (ir->rlistlong <= ir->rvdw))
1223 sprintf(warn_buf, "For energy conservation with switch/shift potentials, %s should be 0.1 to 0.3 nm larger than rvdw.",
1224 IR_TWINRANGE(*ir) ? "rlistlong" : "rlist");
1225 warning_note(wi, warn_buf);
1229 if (ir->vdwtype == evdwUSER && ir->eDispCorr != edispcNO)
1231 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.");
1234 if (ir->nstlist == -1)
1236 sprintf(err_buf, "With nstlist=-1 rvdw and rcoulomb should be smaller than rlist to account for diffusion and possibly charge-group radii");
1237 CHECK(ir->rvdw >= ir->rlist || ir->rcoulomb >= ir->rlist);
1239 sprintf(err_buf, "nstlist can not be smaller than -1");
1240 CHECK(ir->nstlist < -1);
1242 if (ir->eI == eiLBFGS && (ir->coulombtype == eelCUT || ir->vdwtype == evdwCUT)
1245 warning(wi, "For efficient BFGS minimization, use switch/shift/pme instead of cut-off.");
1248 if (ir->eI == eiLBFGS && ir->nbfgscorr <= 0)
1250 warning(wi, "Using L-BFGS with nbfgscorr<=0 just gets you steepest descent.");
1253 /* ENERGY CONSERVATION */
1254 if (ir_NVE(ir) && ir->cutoff_scheme == ecutsGROUP)
1256 if (!EVDW_MIGHT_BE_ZERO_AT_CUTOFF(ir->vdwtype) && ir->rvdw > 0 && ir->vdw_modifier == eintmodNONE)
1258 sprintf(warn_buf, "You are using a cut-off for VdW interactions with NVE, for good energy conservation use vdwtype = %s (possibly with DispCorr)",
1259 evdw_names[evdwSHIFT]);
1260 warning_note(wi, warn_buf);
1262 if (!EEL_MIGHT_BE_ZERO_AT_CUTOFF(ir->coulombtype) && ir->rcoulomb > 0 && ir->coulomb_modifier == eintmodNONE)
1264 sprintf(warn_buf, "You are using a cut-off for electrostatics with NVE, for good energy conservation use coulombtype = %s or %s",
1265 eel_names[eelPMESWITCH], eel_names[eelRF_ZERO]);
1266 warning_note(wi, warn_buf);
1270 /* IMPLICIT SOLVENT */
1271 if (ir->coulombtype == eelGB_NOTUSED)
1273 ir->coulombtype = eelCUT;
1274 ir->implicit_solvent = eisGBSA;
1275 fprintf(stderr, "Note: Old option for generalized born electrostatics given:\n"
1276 "Changing coulombtype from \"generalized-born\" to \"cut-off\" and instead\n"
1277 "setting implicit-solvent value to \"GBSA\" in input section.\n");
1280 if (ir->sa_algorithm == esaSTILL)
1282 sprintf(err_buf, "Still SA algorithm not available yet, use %s or %s instead\n", esa_names[esaAPPROX], esa_names[esaNO]);
1283 CHECK(ir->sa_algorithm == esaSTILL);
1286 if (ir->implicit_solvent == eisGBSA)
1288 sprintf(err_buf, "With GBSA implicit solvent, rgbradii must be equal to rlist.");
1289 CHECK(ir->rgbradii != ir->rlist);
1291 if (ir->coulombtype != eelCUT)
1293 sprintf(err_buf, "With GBSA, coulombtype must be equal to %s\n", eel_names[eelCUT]);
1294 CHECK(ir->coulombtype != eelCUT);
1296 if (ir->vdwtype != evdwCUT)
1298 sprintf(err_buf, "With GBSA, vdw-type must be equal to %s\n", evdw_names[evdwCUT]);
1299 CHECK(ir->vdwtype != evdwCUT);
1301 if (ir->nstgbradii < 1)
1303 sprintf(warn_buf, "Using GBSA with nstgbradii<1, setting nstgbradii=1");
1304 warning_note(wi, warn_buf);
1307 if (ir->sa_algorithm == esaNO)
1309 sprintf(warn_buf, "No SA (non-polar) calculation requested together with GB. Are you sure this is what you want?\n");
1310 warning_note(wi, warn_buf);
1312 if (ir->sa_surface_tension < 0 && ir->sa_algorithm != esaNO)
1314 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");
1315 warning_note(wi, warn_buf);
1317 if (ir->gb_algorithm == egbSTILL)
1319 ir->sa_surface_tension = 0.0049 * CAL2JOULE * 100;
1323 ir->sa_surface_tension = 0.0054 * CAL2JOULE * 100;
1326 if (ir->sa_surface_tension == 0 && ir->sa_algorithm != esaNO)
1328 sprintf(err_buf, "Surface tension set to 0 while SA-calculation requested\n");
1329 CHECK(ir->sa_surface_tension == 0 && ir->sa_algorithm != esaNO);
1336 if (ir->cutoff_scheme != ecutsGROUP)
1338 warning_error(wi, "AdresS simulation supports only cutoff-scheme=group");
1342 warning_error(wi, "AdresS simulation supports only stochastic dynamics");
1344 if (ir->epc != epcNO)
1346 warning_error(wi, "AdresS simulation does not support pressure coupling");
1348 if (EEL_FULL(ir->coulombtype))
1350 warning_error(wi, "AdresS simulation does not support long-range electrostatics");
1355 /* count the number of text elemets separated by whitespace in a string.
1356 str = the input string
1357 maxptr = the maximum number of allowed elements
1358 ptr = the output array of pointers to the first character of each element
1359 returns: the number of elements. */
1360 int str_nelem(const char *str, int maxptr, char *ptr[])
1365 copy0 = strdup(str);
1368 while (*copy != '\0')
1372 gmx_fatal(FARGS, "Too many groups on line: '%s' (max is %d)",
1380 while ((*copy != '\0') && !isspace(*copy))
1399 /* interpret a number of doubles from a string and put them in an array,
1400 after allocating space for them.
1401 str = the input string
1402 n = the (pre-allocated) number of doubles read
1403 r = the output array of doubles. */
1404 static void parse_n_real(char *str, int *n, real **r)
1409 *n = str_nelem(str, MAXPTR, ptr);
1412 for (i = 0; i < *n; i++)
1414 (*r)[i] = strtod(ptr[i], NULL);
1418 static void do_fep_params(t_inputrec *ir, char fep_lambda[][STRLEN], char weights[STRLEN])
1421 int i, j, max_n_lambda, nweights, nfep[efptNR];
1422 t_lambda *fep = ir->fepvals;
1423 t_expanded *expand = ir->expandedvals;
1424 real **count_fep_lambdas;
1425 gmx_bool bOneLambda = TRUE;
1427 snew(count_fep_lambdas, efptNR);
1429 /* FEP input processing */
1430 /* first, identify the number of lambda values for each type.
1431 All that are nonzero must have the same number */
1433 for (i = 0; i < efptNR; i++)
1435 parse_n_real(fep_lambda[i], &(nfep[i]), &(count_fep_lambdas[i]));
1438 /* now, determine the number of components. All must be either zero, or equal. */
1441 for (i = 0; i < efptNR; i++)
1443 if (nfep[i] > max_n_lambda)
1445 max_n_lambda = nfep[i]; /* here's a nonzero one. All of them
1446 must have the same number if its not zero.*/
1451 for (i = 0; i < efptNR; i++)
1455 ir->fepvals->separate_dvdl[i] = FALSE;
1457 else if (nfep[i] == max_n_lambda)
1459 if (i != efptTEMPERATURE) /* we treat this differently -- not really a reason to compute the derivative with
1460 respect to the temperature currently */
1462 ir->fepvals->separate_dvdl[i] = TRUE;
1467 gmx_fatal(FARGS, "Number of lambdas (%d) for FEP type %s not equal to number of other types (%d)",
1468 nfep[i], efpt_names[i], max_n_lambda);
1471 /* we don't print out dhdl if the temperature is changing, since we can't correctly define dhdl in this case */
1472 ir->fepvals->separate_dvdl[efptTEMPERATURE] = FALSE;
1474 /* the number of lambdas is the number we've read in, which is either zero
1475 or the same for all */
1476 fep->n_lambda = max_n_lambda;
1478 /* allocate space for the array of lambda values */
1479 snew(fep->all_lambda, efptNR);
1480 /* if init_lambda is defined, we need to set lambda */
1481 if ((fep->init_lambda > 0) && (fep->n_lambda == 0))
1483 ir->fepvals->separate_dvdl[efptFEP] = TRUE;
1485 /* otherwise allocate the space for all of the lambdas, and transfer the data */
1486 for (i = 0; i < efptNR; i++)
1488 snew(fep->all_lambda[i], fep->n_lambda);
1489 if (nfep[i] > 0) /* if it's zero, then the count_fep_lambda arrays
1492 for (j = 0; j < fep->n_lambda; j++)
1494 fep->all_lambda[i][j] = (double)count_fep_lambdas[i][j];
1496 sfree(count_fep_lambdas[i]);
1499 sfree(count_fep_lambdas);
1501 /* "fep-vals" is either zero or the full number. If zero, we'll need to define fep-lambdas for internal
1502 bookkeeping -- for now, init_lambda */
1504 if ((nfep[efptFEP] == 0) && (fep->init_lambda >= 0))
1506 for (i = 0; i < fep->n_lambda; i++)
1508 fep->all_lambda[efptFEP][i] = fep->init_lambda;
1512 /* check to see if only a single component lambda is defined, and soft core is defined.
1513 In this case, turn on coulomb soft core */
1515 if (max_n_lambda == 0)
1521 for (i = 0; i < efptNR; i++)
1523 if ((nfep[i] != 0) && (i != efptFEP))
1529 if ((bOneLambda) && (fep->sc_alpha > 0))
1531 fep->bScCoul = TRUE;
1534 /* Fill in the others with the efptFEP if they are not explicitly
1535 specified (i.e. nfep[i] == 0). This means if fep is not defined,
1536 they are all zero. */
1538 for (i = 0; i < efptNR; i++)
1540 if ((nfep[i] == 0) && (i != efptFEP))
1542 for (j = 0; j < fep->n_lambda; j++)
1544 fep->all_lambda[i][j] = fep->all_lambda[efptFEP][j];
1550 /* make it easier if sc_r_power = 48 by increasing it to the 4th power, to be in the right scale. */
1551 if (fep->sc_r_power == 48)
1553 if (fep->sc_alpha > 0.1)
1555 gmx_fatal(FARGS, "sc_alpha (%f) for sc_r_power = 48 should usually be between 0.001 and 0.004", fep->sc_alpha);
1559 expand = ir->expandedvals;
1560 /* now read in the weights */
1561 parse_n_real(weights, &nweights, &(expand->init_lambda_weights));
1564 snew(expand->init_lambda_weights, fep->n_lambda); /* initialize to zero */
1566 else if (nweights != fep->n_lambda)
1568 gmx_fatal(FARGS, "Number of weights (%d) is not equal to number of lambda values (%d)",
1569 nweights, fep->n_lambda);
1571 if ((expand->nstexpanded < 0) && (ir->efep != efepNO))
1573 expand->nstexpanded = fep->nstdhdl;
1574 /* if you don't specify nstexpanded when doing expanded ensemble free energy calcs, it is set to nstdhdl */
1576 if ((expand->nstexpanded < 0) && ir->bSimTemp)
1578 expand->nstexpanded = 2*(int)(ir->opts.tau_t[0]/ir->delta_t);
1579 /* if you don't specify nstexpanded when doing expanded ensemble simulated tempering, it is set to
1580 2*tau_t just to be careful so it's not to frequent */
1585 static void do_simtemp_params(t_inputrec *ir)
1588 snew(ir->simtempvals->temperatures, ir->fepvals->n_lambda);
1589 GetSimTemps(ir->fepvals->n_lambda, ir->simtempvals, ir->fepvals->all_lambda[efptTEMPERATURE]);
1594 static void do_wall_params(t_inputrec *ir,
1595 char *wall_atomtype, char *wall_density,
1599 char *names[MAXPTR];
1602 opts->wall_atomtype[0] = NULL;
1603 opts->wall_atomtype[1] = NULL;
1605 ir->wall_atomtype[0] = -1;
1606 ir->wall_atomtype[1] = -1;
1607 ir->wall_density[0] = 0;
1608 ir->wall_density[1] = 0;
1612 nstr = str_nelem(wall_atomtype, MAXPTR, names);
1613 if (nstr != ir->nwall)
1615 gmx_fatal(FARGS, "Expected %d elements for wall_atomtype, found %d",
1618 for (i = 0; i < ir->nwall; i++)
1620 opts->wall_atomtype[i] = strdup(names[i]);
1623 if (ir->wall_type == ewt93 || ir->wall_type == ewt104)
1625 nstr = str_nelem(wall_density, MAXPTR, names);
1626 if (nstr != ir->nwall)
1628 gmx_fatal(FARGS, "Expected %d elements for wall-density, found %d", ir->nwall, nstr);
1630 for (i = 0; i < ir->nwall; i++)
1632 sscanf(names[i], "%lf", &dbl);
1635 gmx_fatal(FARGS, "wall-density[%d] = %f\n", i, dbl);
1637 ir->wall_density[i] = dbl;
1643 static void add_wall_energrps(gmx_groups_t *groups, int nwall, t_symtab *symtab)
1651 srenew(groups->grpname, groups->ngrpname+nwall);
1652 grps = &(groups->grps[egcENER]);
1653 srenew(grps->nm_ind, grps->nr+nwall);
1654 for (i = 0; i < nwall; i++)
1656 sprintf(str, "wall%d", i);
1657 groups->grpname[groups->ngrpname] = put_symtab(symtab, str);
1658 grps->nm_ind[grps->nr++] = groups->ngrpname++;
1663 void read_expandedparams(int *ninp_p, t_inpfile **inp_p,
1664 t_expanded *expand, warninp_t wi)
1666 int ninp, nerror = 0;
1672 /* read expanded ensemble parameters */
1673 CCTYPE ("expanded ensemble variables");
1674 ITYPE ("nstexpanded", expand->nstexpanded, -1);
1675 EETYPE("lmc-stats", expand->elamstats, elamstats_names);
1676 EETYPE("lmc-move", expand->elmcmove, elmcmove_names);
1677 EETYPE("lmc-weights-equil", expand->elmceq, elmceq_names);
1678 ITYPE ("weight-equil-number-all-lambda", expand->equil_n_at_lam, -1);
1679 ITYPE ("weight-equil-number-samples", expand->equil_samples, -1);
1680 ITYPE ("weight-equil-number-steps", expand->equil_steps, -1);
1681 RTYPE ("weight-equil-wl-delta", expand->equil_wl_delta, -1);
1682 RTYPE ("weight-equil-count-ratio", expand->equil_ratio, -1);
1683 CCTYPE("Seed for Monte Carlo in lambda space");
1684 ITYPE ("lmc-seed", expand->lmc_seed, -1);
1685 RTYPE ("mc-temperature", expand->mc_temp, -1);
1686 ITYPE ("lmc-repeats", expand->lmc_repeats, 1);
1687 ITYPE ("lmc-gibbsdelta", expand->gibbsdeltalam, -1);
1688 ITYPE ("lmc-forced-nstart", expand->lmc_forced_nstart, 0);
1689 EETYPE("symmetrized-transition-matrix", expand->bSymmetrizedTMatrix, yesno_names);
1690 ITYPE("nst-transition-matrix", expand->nstTij, -1);
1691 ITYPE ("mininum-var-min", expand->minvarmin, 100); /*default is reasonable */
1692 ITYPE ("weight-c-range", expand->c_range, 0); /* default is just C=0 */
1693 RTYPE ("wl-scale", expand->wl_scale, 0.8);
1694 RTYPE ("wl-ratio", expand->wl_ratio, 0.8);
1695 RTYPE ("init-wl-delta", expand->init_wl_delta, 1.0);
1696 EETYPE("wl-oneovert", expand->bWLoneovert, yesno_names);
1704 void get_ir(const char *mdparin, const char *mdparout,
1705 t_inputrec *ir, t_gromppopts *opts,
1709 double dumdub[2][6];
1713 char warn_buf[STRLEN];
1714 t_lambda *fep = ir->fepvals;
1715 t_expanded *expand = ir->expandedvals;
1717 init_inputrec_strings();
1718 inp = read_inpfile(mdparin, &ninp, wi);
1720 snew(dumstr[0], STRLEN);
1721 snew(dumstr[1], STRLEN);
1723 if (-1 == search_einp(ninp, inp, "cutoff-scheme"))
1726 "%s did not specify a value for the .mdp option "
1727 "\"cutoff-scheme\". Probably it was first intended for use "
1728 "with GROMACS before 4.6. In 4.6, the Verlet scheme was "
1729 "introduced, but the group scheme was still the default. "
1730 "The default is now the Verlet scheme, so you will observe "
1731 "different behaviour.", mdparin);
1732 warning_note(wi, warn_buf);
1735 /* remove the following deprecated commands */
1738 REM_TYPE("domain-decomposition");
1739 REM_TYPE("andersen-seed");
1741 REM_TYPE("dihre-fc");
1742 REM_TYPE("dihre-tau");
1743 REM_TYPE("nstdihreout");
1744 REM_TYPE("nstcheckpoint");
1746 /* replace the following commands with the clearer new versions*/
1747 REPL_TYPE("unconstrained-start", "continuation");
1748 REPL_TYPE("foreign-lambda", "fep-lambdas");
1749 REPL_TYPE("verlet-buffer-drift", "verlet-buffer-tolerance");
1751 CCTYPE ("VARIOUS PREPROCESSING OPTIONS");
1752 CTYPE ("Preprocessor information: use cpp syntax.");
1753 CTYPE ("e.g.: -I/home/joe/doe -I/home/mary/roe");
1754 STYPE ("include", opts->include, NULL);
1755 CTYPE ("e.g.: -DPOSRES -DFLEXIBLE (note these variable names are case sensitive)");
1756 STYPE ("define", opts->define, NULL);
1758 CCTYPE ("RUN CONTROL PARAMETERS");
1759 EETYPE("integrator", ir->eI, ei_names);
1760 CTYPE ("Start time and timestep in ps");
1761 RTYPE ("tinit", ir->init_t, 0.0);
1762 RTYPE ("dt", ir->delta_t, 0.001);
1763 STEPTYPE ("nsteps", ir->nsteps, 0);
1764 CTYPE ("For exact run continuation or redoing part of a run");
1765 STEPTYPE ("init-step", ir->init_step, 0);
1766 CTYPE ("Part index is updated automatically on checkpointing (keeps files separate)");
1767 ITYPE ("simulation-part", ir->simulation_part, 1);
1768 CTYPE ("mode for center of mass motion removal");
1769 EETYPE("comm-mode", ir->comm_mode, ecm_names);
1770 CTYPE ("number of steps for center of mass motion removal");
1771 ITYPE ("nstcomm", ir->nstcomm, 100);
1772 CTYPE ("group(s) for center of mass motion removal");
1773 STYPE ("comm-grps", is->vcm, NULL);
1775 CCTYPE ("LANGEVIN DYNAMICS OPTIONS");
1776 CTYPE ("Friction coefficient (amu/ps) and random seed");
1777 RTYPE ("bd-fric", ir->bd_fric, 0.0);
1778 ITYPE ("ld-seed", ir->ld_seed, 1993);
1781 CCTYPE ("ENERGY MINIMIZATION OPTIONS");
1782 CTYPE ("Force tolerance and initial step-size");
1783 RTYPE ("emtol", ir->em_tol, 10.0);
1784 RTYPE ("emstep", ir->em_stepsize, 0.01);
1785 CTYPE ("Max number of iterations in relax-shells");
1786 ITYPE ("niter", ir->niter, 20);
1787 CTYPE ("Step size (ps^2) for minimization of flexible constraints");
1788 RTYPE ("fcstep", ir->fc_stepsize, 0);
1789 CTYPE ("Frequency of steepest descents steps when doing CG");
1790 ITYPE ("nstcgsteep", ir->nstcgsteep, 1000);
1791 ITYPE ("nbfgscorr", ir->nbfgscorr, 10);
1793 CCTYPE ("TEST PARTICLE INSERTION OPTIONS");
1794 RTYPE ("rtpi", ir->rtpi, 0.05);
1796 /* Output options */
1797 CCTYPE ("OUTPUT CONTROL OPTIONS");
1798 CTYPE ("Output frequency for coords (x), velocities (v) and forces (f)");
1799 ITYPE ("nstxout", ir->nstxout, 0);
1800 ITYPE ("nstvout", ir->nstvout, 0);
1801 ITYPE ("nstfout", ir->nstfout, 0);
1802 ir->nstcheckpoint = 1000;
1803 CTYPE ("Output frequency for energies to log file and energy file");
1804 ITYPE ("nstlog", ir->nstlog, 1000);
1805 ITYPE ("nstcalcenergy", ir->nstcalcenergy, 100);
1806 ITYPE ("nstenergy", ir->nstenergy, 1000);
1807 CTYPE ("Output frequency and precision for .xtc file");
1808 ITYPE ("nstxtcout", ir->nstxtcout, 0);
1809 RTYPE ("xtc-precision", ir->xtcprec, 1000.0);
1810 CTYPE ("This selects the subset of atoms for the .xtc file. You can");
1811 CTYPE ("select multiple groups. By default all atoms will be written.");
1812 STYPE ("xtc-grps", is->xtc_grps, NULL);
1813 CTYPE ("Selection of energy groups");
1814 STYPE ("energygrps", is->energy, NULL);
1816 /* Neighbor searching */
1817 CCTYPE ("NEIGHBORSEARCHING PARAMETERS");
1818 CTYPE ("cut-off scheme (Verlet: particle based cut-offs, group: using charge groups)");
1819 EETYPE("cutoff-scheme", ir->cutoff_scheme, ecutscheme_names);
1820 CTYPE ("nblist update frequency");
1821 ITYPE ("nstlist", ir->nstlist, 10);
1822 CTYPE ("ns algorithm (simple or grid)");
1823 EETYPE("ns-type", ir->ns_type, ens_names);
1824 /* set ndelta to the optimal value of 2 */
1826 CTYPE ("Periodic boundary conditions: xyz, no, xy");
1827 EETYPE("pbc", ir->ePBC, epbc_names);
1828 EETYPE("periodic-molecules", ir->bPeriodicMols, yesno_names);
1829 CTYPE ("Allowed energy error due to the Verlet buffer in kJ/mol/ps per atom,");
1830 CTYPE ("a value of -1 means: use rlist");
1831 RTYPE("verlet-buffer-tolerance", ir->verletbuf_tol, 0.005);
1832 CTYPE ("nblist cut-off");
1833 RTYPE ("rlist", ir->rlist, 1.0);
1834 CTYPE ("long-range cut-off for switched potentials");
1835 RTYPE ("rlistlong", ir->rlistlong, -1);
1836 ITYPE ("nstcalclr", ir->nstcalclr, -1);
1838 /* Electrostatics */
1839 CCTYPE ("OPTIONS FOR ELECTROSTATICS AND VDW");
1840 CTYPE ("Method for doing electrostatics");
1841 EETYPE("coulombtype", ir->coulombtype, eel_names);
1842 EETYPE("coulomb-modifier", ir->coulomb_modifier, eintmod_names);
1843 CTYPE ("cut-off lengths");
1844 RTYPE ("rcoulomb-switch", ir->rcoulomb_switch, 0.0);
1845 RTYPE ("rcoulomb", ir->rcoulomb, 1.0);
1846 CTYPE ("Relative dielectric constant for the medium and the reaction field");
1847 RTYPE ("epsilon-r", ir->epsilon_r, 1.0);
1848 RTYPE ("epsilon-rf", ir->epsilon_rf, 0.0);
1849 CTYPE ("Method for doing Van der Waals");
1850 EETYPE("vdw-type", ir->vdwtype, evdw_names);
1851 EETYPE("vdw-modifier", ir->vdw_modifier, eintmod_names);
1852 CTYPE ("cut-off lengths");
1853 RTYPE ("rvdw-switch", ir->rvdw_switch, 0.0);
1854 RTYPE ("rvdw", ir->rvdw, 1.0);
1855 CTYPE ("Apply long range dispersion corrections for Energy and Pressure");
1856 EETYPE("DispCorr", ir->eDispCorr, edispc_names);
1857 CTYPE ("Extension of the potential lookup tables beyond the cut-off");
1858 RTYPE ("table-extension", ir->tabext, 1.0);
1859 CTYPE ("Separate tables between energy group pairs");
1860 STYPE ("energygrp-table", is->egptable, NULL);
1861 CTYPE ("Spacing for the PME/PPPM FFT grid");
1862 RTYPE ("fourierspacing", ir->fourier_spacing, 0.12);
1863 CTYPE ("FFT grid size, when a value is 0 fourierspacing will be used");
1864 ITYPE ("fourier-nx", ir->nkx, 0);
1865 ITYPE ("fourier-ny", ir->nky, 0);
1866 ITYPE ("fourier-nz", ir->nkz, 0);
1867 CTYPE ("EWALD/PME/PPPM parameters");
1868 ITYPE ("pme-order", ir->pme_order, 4);
1869 RTYPE ("ewald-rtol", ir->ewald_rtol, 0.00001);
1870 RTYPE ("ewald-rtol-lj", ir->ewald_rtol_lj, 0.001);
1871 EETYPE("lj-pme-comb-rule", ir->ljpme_combination_rule, eljpme_names);
1872 EETYPE("ewald-geometry", ir->ewald_geometry, eewg_names);
1873 RTYPE ("epsilon-surface", ir->epsilon_surface, 0.0);
1874 EETYPE("optimize-fft", ir->bOptFFT, yesno_names);
1876 CCTYPE("IMPLICIT SOLVENT ALGORITHM");
1877 EETYPE("implicit-solvent", ir->implicit_solvent, eis_names);
1879 CCTYPE ("GENERALIZED BORN ELECTROSTATICS");
1880 CTYPE ("Algorithm for calculating Born radii");
1881 EETYPE("gb-algorithm", ir->gb_algorithm, egb_names);
1882 CTYPE ("Frequency of calculating the Born radii inside rlist");
1883 ITYPE ("nstgbradii", ir->nstgbradii, 1);
1884 CTYPE ("Cutoff for Born radii calculation; the contribution from atoms");
1885 CTYPE ("between rlist and rgbradii is updated every nstlist steps");
1886 RTYPE ("rgbradii", ir->rgbradii, 1.0);
1887 CTYPE ("Dielectric coefficient of the implicit solvent");
1888 RTYPE ("gb-epsilon-solvent", ir->gb_epsilon_solvent, 80.0);
1889 CTYPE ("Salt concentration in M for Generalized Born models");
1890 RTYPE ("gb-saltconc", ir->gb_saltconc, 0.0);
1891 CTYPE ("Scaling factors used in the OBC GB model. Default values are OBC(II)");
1892 RTYPE ("gb-obc-alpha", ir->gb_obc_alpha, 1.0);
1893 RTYPE ("gb-obc-beta", ir->gb_obc_beta, 0.8);
1894 RTYPE ("gb-obc-gamma", ir->gb_obc_gamma, 4.85);
1895 RTYPE ("gb-dielectric-offset", ir->gb_dielectric_offset, 0.009);
1896 EETYPE("sa-algorithm", ir->sa_algorithm, esa_names);
1897 CTYPE ("Surface tension (kJ/mol/nm^2) for the SA (nonpolar surface) part of GBSA");
1898 CTYPE ("The value -1 will set default value for Still/HCT/OBC GB-models.");
1899 RTYPE ("sa-surface-tension", ir->sa_surface_tension, -1);
1901 /* Coupling stuff */
1902 CCTYPE ("OPTIONS FOR WEAK COUPLING ALGORITHMS");
1903 CTYPE ("Temperature coupling");
1904 EETYPE("tcoupl", ir->etc, etcoupl_names);
1905 ITYPE ("nsttcouple", ir->nsttcouple, -1);
1906 ITYPE("nh-chain-length", ir->opts.nhchainlength, 10);
1907 EETYPE("print-nose-hoover-chain-variables", ir->bPrintNHChains, yesno_names);
1908 CTYPE ("Groups to couple separately");
1909 STYPE ("tc-grps", is->tcgrps, NULL);
1910 CTYPE ("Time constant (ps) and reference temperature (K)");
1911 STYPE ("tau-t", is->tau_t, NULL);
1912 STYPE ("ref-t", is->ref_t, NULL);
1913 CTYPE ("pressure coupling");
1914 EETYPE("pcoupl", ir->epc, epcoupl_names);
1915 EETYPE("pcoupltype", ir->epct, epcoupltype_names);
1916 ITYPE ("nstpcouple", ir->nstpcouple, -1);
1917 CTYPE ("Time constant (ps), compressibility (1/bar) and reference P (bar)");
1918 RTYPE ("tau-p", ir->tau_p, 1.0);
1919 STYPE ("compressibility", dumstr[0], NULL);
1920 STYPE ("ref-p", dumstr[1], NULL);
1921 CTYPE ("Scaling of reference coordinates, No, All or COM");
1922 EETYPE ("refcoord-scaling", ir->refcoord_scaling, erefscaling_names);
1925 CCTYPE ("OPTIONS FOR QMMM calculations");
1926 EETYPE("QMMM", ir->bQMMM, yesno_names);
1927 CTYPE ("Groups treated Quantum Mechanically");
1928 STYPE ("QMMM-grps", is->QMMM, NULL);
1929 CTYPE ("QM method");
1930 STYPE("QMmethod", is->QMmethod, NULL);
1931 CTYPE ("QMMM scheme");
1932 EETYPE("QMMMscheme", ir->QMMMscheme, eQMMMscheme_names);
1933 CTYPE ("QM basisset");
1934 STYPE("QMbasis", is->QMbasis, NULL);
1935 CTYPE ("QM charge");
1936 STYPE ("QMcharge", is->QMcharge, NULL);
1937 CTYPE ("QM multiplicity");
1938 STYPE ("QMmult", is->QMmult, NULL);
1939 CTYPE ("Surface Hopping");
1940 STYPE ("SH", is->bSH, NULL);
1941 CTYPE ("CAS space options");
1942 STYPE ("CASorbitals", is->CASorbitals, NULL);
1943 STYPE ("CASelectrons", is->CASelectrons, NULL);
1944 STYPE ("SAon", is->SAon, NULL);
1945 STYPE ("SAoff", is->SAoff, NULL);
1946 STYPE ("SAsteps", is->SAsteps, NULL);
1947 CTYPE ("Scale factor for MM charges");
1948 RTYPE ("MMChargeScaleFactor", ir->scalefactor, 1.0);
1949 CTYPE ("Optimization of QM subsystem");
1950 STYPE ("bOPT", is->bOPT, NULL);
1951 STYPE ("bTS", is->bTS, NULL);
1953 /* Simulated annealing */
1954 CCTYPE("SIMULATED ANNEALING");
1955 CTYPE ("Type of annealing for each temperature group (no/single/periodic)");
1956 STYPE ("annealing", is->anneal, NULL);
1957 CTYPE ("Number of time points to use for specifying annealing in each group");
1958 STYPE ("annealing-npoints", is->anneal_npoints, NULL);
1959 CTYPE ("List of times at the annealing points for each group");
1960 STYPE ("annealing-time", is->anneal_time, NULL);
1961 CTYPE ("Temp. at each annealing point, for each group.");
1962 STYPE ("annealing-temp", is->anneal_temp, NULL);
1965 CCTYPE ("GENERATE VELOCITIES FOR STARTUP RUN");
1966 EETYPE("gen-vel", opts->bGenVel, yesno_names);
1967 RTYPE ("gen-temp", opts->tempi, 300.0);
1968 ITYPE ("gen-seed", opts->seed, 173529);
1971 CCTYPE ("OPTIONS FOR BONDS");
1972 EETYPE("constraints", opts->nshake, constraints);
1973 CTYPE ("Type of constraint algorithm");
1974 EETYPE("constraint-algorithm", ir->eConstrAlg, econstr_names);
1975 CTYPE ("Do not constrain the start configuration");
1976 EETYPE("continuation", ir->bContinuation, yesno_names);
1977 CTYPE ("Use successive overrelaxation to reduce the number of shake iterations");
1978 EETYPE("Shake-SOR", ir->bShakeSOR, yesno_names);
1979 CTYPE ("Relative tolerance of shake");
1980 RTYPE ("shake-tol", ir->shake_tol, 0.0001);
1981 CTYPE ("Highest order in the expansion of the constraint coupling matrix");
1982 ITYPE ("lincs-order", ir->nProjOrder, 4);
1983 CTYPE ("Number of iterations in the final step of LINCS. 1 is fine for");
1984 CTYPE ("normal simulations, but use 2 to conserve energy in NVE runs.");
1985 CTYPE ("For energy minimization with constraints it should be 4 to 8.");
1986 ITYPE ("lincs-iter", ir->nLincsIter, 1);
1987 CTYPE ("Lincs will write a warning to the stderr if in one step a bond");
1988 CTYPE ("rotates over more degrees than");
1989 RTYPE ("lincs-warnangle", ir->LincsWarnAngle, 30.0);
1990 CTYPE ("Convert harmonic bonds to morse potentials");
1991 EETYPE("morse", opts->bMorse, yesno_names);
1993 /* Energy group exclusions */
1994 CCTYPE ("ENERGY GROUP EXCLUSIONS");
1995 CTYPE ("Pairs of energy groups for which all non-bonded interactions are excluded");
1996 STYPE ("energygrp-excl", is->egpexcl, NULL);
2000 CTYPE ("Number of walls, type, atom types, densities and box-z scale factor for Ewald");
2001 ITYPE ("nwall", ir->nwall, 0);
2002 EETYPE("wall-type", ir->wall_type, ewt_names);
2003 RTYPE ("wall-r-linpot", ir->wall_r_linpot, -1);
2004 STYPE ("wall-atomtype", is->wall_atomtype, NULL);
2005 STYPE ("wall-density", is->wall_density, NULL);
2006 RTYPE ("wall-ewald-zfac", ir->wall_ewald_zfac, 3);
2009 CCTYPE("COM PULLING");
2010 CTYPE("Pull type: no, umbrella, constraint or constant-force");
2011 EETYPE("pull", ir->ePull, epull_names);
2012 if (ir->ePull != epullNO)
2015 is->pull_grp = read_pullparams(&ninp, &inp, ir->pull, &opts->pull_start, wi);
2018 /* Enforced rotation */
2019 CCTYPE("ENFORCED ROTATION");
2020 CTYPE("Enforced rotation: No or Yes");
2021 EETYPE("rotation", ir->bRot, yesno_names);
2025 is->rot_grp = read_rotparams(&ninp, &inp, ir->rot, wi);
2029 CCTYPE("NMR refinement stuff");
2030 CTYPE ("Distance restraints type: No, Simple or Ensemble");
2031 EETYPE("disre", ir->eDisre, edisre_names);
2032 CTYPE ("Force weighting of pairs in one distance restraint: Conservative or Equal");
2033 EETYPE("disre-weighting", ir->eDisreWeighting, edisreweighting_names);
2034 CTYPE ("Use sqrt of the time averaged times the instantaneous violation");
2035 EETYPE("disre-mixed", ir->bDisreMixed, yesno_names);
2036 RTYPE ("disre-fc", ir->dr_fc, 1000.0);
2037 RTYPE ("disre-tau", ir->dr_tau, 0.0);
2038 CTYPE ("Output frequency for pair distances to energy file");
2039 ITYPE ("nstdisreout", ir->nstdisreout, 100);
2040 CTYPE ("Orientation restraints: No or Yes");
2041 EETYPE("orire", opts->bOrire, yesno_names);
2042 CTYPE ("Orientation restraints force constant and tau for time averaging");
2043 RTYPE ("orire-fc", ir->orires_fc, 0.0);
2044 RTYPE ("orire-tau", ir->orires_tau, 0.0);
2045 STYPE ("orire-fitgrp", is->orirefitgrp, NULL);
2046 CTYPE ("Output frequency for trace(SD) and S to energy file");
2047 ITYPE ("nstorireout", ir->nstorireout, 100);
2049 /* free energy variables */
2050 CCTYPE ("Free energy variables");
2051 EETYPE("free-energy", ir->efep, efep_names);
2052 STYPE ("couple-moltype", is->couple_moltype, NULL);
2053 EETYPE("couple-lambda0", opts->couple_lam0, couple_lam);
2054 EETYPE("couple-lambda1", opts->couple_lam1, couple_lam);
2055 EETYPE("couple-intramol", opts->bCoupleIntra, yesno_names);
2057 RTYPE ("init-lambda", fep->init_lambda, -1); /* start with -1 so
2059 it was not entered */
2060 ITYPE ("init-lambda-state", fep->init_fep_state, -1);
2061 RTYPE ("delta-lambda", fep->delta_lambda, 0.0);
2062 ITYPE ("nstdhdl", fep->nstdhdl, 50);
2063 STYPE ("fep-lambdas", is->fep_lambda[efptFEP], NULL);
2064 STYPE ("mass-lambdas", is->fep_lambda[efptMASS], NULL);
2065 STYPE ("coul-lambdas", is->fep_lambda[efptCOUL], NULL);
2066 STYPE ("vdw-lambdas", is->fep_lambda[efptVDW], NULL);
2067 STYPE ("bonded-lambdas", is->fep_lambda[efptBONDED], NULL);
2068 STYPE ("restraint-lambdas", is->fep_lambda[efptRESTRAINT], NULL);
2069 STYPE ("temperature-lambdas", is->fep_lambda[efptTEMPERATURE], NULL);
2070 ITYPE ("calc-lambda-neighbors", fep->lambda_neighbors, 1);
2071 STYPE ("init-lambda-weights", is->lambda_weights, NULL);
2072 EETYPE("dhdl-print-energy", fep->bPrintEnergy, yesno_names);
2073 RTYPE ("sc-alpha", fep->sc_alpha, 0.0);
2074 ITYPE ("sc-power", fep->sc_power, 1);
2075 RTYPE ("sc-r-power", fep->sc_r_power, 6.0);
2076 RTYPE ("sc-sigma", fep->sc_sigma, 0.3);
2077 EETYPE("sc-coul", fep->bScCoul, yesno_names);
2078 ITYPE ("dh_hist_size", fep->dh_hist_size, 0);
2079 RTYPE ("dh_hist_spacing", fep->dh_hist_spacing, 0.1);
2080 EETYPE("separate-dhdl-file", fep->separate_dhdl_file,
2081 separate_dhdl_file_names);
2082 EETYPE("dhdl-derivatives", fep->dhdl_derivatives, dhdl_derivatives_names);
2083 ITYPE ("dh_hist_size", fep->dh_hist_size, 0);
2084 RTYPE ("dh_hist_spacing", fep->dh_hist_spacing, 0.1);
2086 /* Non-equilibrium MD stuff */
2087 CCTYPE("Non-equilibrium MD stuff");
2088 STYPE ("acc-grps", is->accgrps, NULL);
2089 STYPE ("accelerate", is->acc, NULL);
2090 STYPE ("freezegrps", is->freeze, NULL);
2091 STYPE ("freezedim", is->frdim, NULL);
2092 RTYPE ("cos-acceleration", ir->cos_accel, 0);
2093 STYPE ("deform", is->deform, NULL);
2095 /* simulated tempering variables */
2096 CCTYPE("simulated tempering variables");
2097 EETYPE("simulated-tempering", ir->bSimTemp, yesno_names);
2098 EETYPE("simulated-tempering-scaling", ir->simtempvals->eSimTempScale, esimtemp_names);
2099 RTYPE("sim-temp-low", ir->simtempvals->simtemp_low, 300.0);
2100 RTYPE("sim-temp-high", ir->simtempvals->simtemp_high, 300.0);
2102 /* expanded ensemble variables */
2103 if (ir->efep == efepEXPANDED || ir->bSimTemp)
2105 read_expandedparams(&ninp, &inp, expand, wi);
2108 /* Electric fields */
2109 CCTYPE("Electric fields");
2110 CTYPE ("Format is number of terms (int) and for all terms an amplitude (real)");
2111 CTYPE ("and a phase angle (real)");
2112 STYPE ("E-x", is->efield_x, NULL);
2113 STYPE ("E-xt", is->efield_xt, NULL);
2114 STYPE ("E-y", is->efield_y, NULL);
2115 STYPE ("E-yt", is->efield_yt, NULL);
2116 STYPE ("E-z", is->efield_z, NULL);
2117 STYPE ("E-zt", is->efield_zt, NULL);
2119 /* AdResS defined thingies */
2120 CCTYPE ("AdResS parameters");
2121 EETYPE("adress", ir->bAdress, yesno_names);
2124 snew(ir->adress, 1);
2125 read_adressparams(&ninp, &inp, ir->adress, wi);
2128 /* User defined thingies */
2129 CCTYPE ("User defined thingies");
2130 STYPE ("user1-grps", is->user1, NULL);
2131 STYPE ("user2-grps", is->user2, NULL);
2132 ITYPE ("userint1", ir->userint1, 0);
2133 ITYPE ("userint2", ir->userint2, 0);
2134 ITYPE ("userint3", ir->userint3, 0);
2135 ITYPE ("userint4", ir->userint4, 0);
2136 RTYPE ("userreal1", ir->userreal1, 0);
2137 RTYPE ("userreal2", ir->userreal2, 0);
2138 RTYPE ("userreal3", ir->userreal3, 0);
2139 RTYPE ("userreal4", ir->userreal4, 0);
2142 write_inpfile(mdparout, ninp, inp, FALSE, wi);
2143 for (i = 0; (i < ninp); i++)
2146 sfree(inp[i].value);
2150 /* Process options if necessary */
2151 for (m = 0; m < 2; m++)
2153 for (i = 0; i < 2*DIM; i++)
2162 if (sscanf(dumstr[m], "%lf", &(dumdub[m][XX])) != 1)
2164 warning_error(wi, "Pressure coupling not enough values (I need 1)");
2166 dumdub[m][YY] = dumdub[m][ZZ] = dumdub[m][XX];
2168 case epctSEMIISOTROPIC:
2169 case epctSURFACETENSION:
2170 if (sscanf(dumstr[m], "%lf%lf",
2171 &(dumdub[m][XX]), &(dumdub[m][ZZ])) != 2)
2173 warning_error(wi, "Pressure coupling not enough values (I need 2)");
2175 dumdub[m][YY] = dumdub[m][XX];
2177 case epctANISOTROPIC:
2178 if (sscanf(dumstr[m], "%lf%lf%lf%lf%lf%lf",
2179 &(dumdub[m][XX]), &(dumdub[m][YY]), &(dumdub[m][ZZ]),
2180 &(dumdub[m][3]), &(dumdub[m][4]), &(dumdub[m][5])) != 6)
2182 warning_error(wi, "Pressure coupling not enough values (I need 6)");
2186 gmx_fatal(FARGS, "Pressure coupling type %s not implemented yet",
2187 epcoupltype_names[ir->epct]);
2191 clear_mat(ir->ref_p);
2192 clear_mat(ir->compress);
2193 for (i = 0; i < DIM; i++)
2195 ir->ref_p[i][i] = dumdub[1][i];
2196 ir->compress[i][i] = dumdub[0][i];
2198 if (ir->epct == epctANISOTROPIC)
2200 ir->ref_p[XX][YY] = dumdub[1][3];
2201 ir->ref_p[XX][ZZ] = dumdub[1][4];
2202 ir->ref_p[YY][ZZ] = dumdub[1][5];
2203 if (ir->ref_p[XX][YY] != 0 && ir->ref_p[XX][ZZ] != 0 && ir->ref_p[YY][ZZ] != 0)
2205 warning(wi, "All off-diagonal reference pressures are non-zero. Are you sure you want to apply a threefold shear stress?\n");
2207 ir->compress[XX][YY] = dumdub[0][3];
2208 ir->compress[XX][ZZ] = dumdub[0][4];
2209 ir->compress[YY][ZZ] = dumdub[0][5];
2210 for (i = 0; i < DIM; i++)
2212 for (m = 0; m < i; m++)
2214 ir->ref_p[i][m] = ir->ref_p[m][i];
2215 ir->compress[i][m] = ir->compress[m][i];
2220 if (ir->comm_mode == ecmNO)
2225 opts->couple_moltype = NULL;
2226 if (strlen(is->couple_moltype) > 0)
2228 if (ir->efep != efepNO)
2230 opts->couple_moltype = strdup(is->couple_moltype);
2231 if (opts->couple_lam0 == opts->couple_lam1)
2233 warning(wi, "The lambda=0 and lambda=1 states for coupling are identical");
2235 if (ir->eI == eiMD && (opts->couple_lam0 == ecouplamNONE ||
2236 opts->couple_lam1 == ecouplamNONE))
2238 warning(wi, "For proper sampling of the (nearly) decoupled state, stochastic dynamics should be used");
2243 warning(wi, "Can not couple a molecule with free_energy = no");
2246 /* FREE ENERGY AND EXPANDED ENSEMBLE OPTIONS */
2247 if (ir->efep != efepNO)
2249 if (fep->delta_lambda > 0)
2251 ir->efep = efepSLOWGROWTH;
2257 fep->bPrintEnergy = TRUE;
2258 /* always print out the energy to dhdl if we are doing expanded ensemble, since we need the total energy
2259 if the temperature is changing. */
2262 if ((ir->efep != efepNO) || ir->bSimTemp)
2264 ir->bExpanded = FALSE;
2265 if ((ir->efep == efepEXPANDED) || ir->bSimTemp)
2267 ir->bExpanded = TRUE;
2269 do_fep_params(ir, is->fep_lambda, is->lambda_weights);
2270 if (ir->bSimTemp) /* done after fep params */
2272 do_simtemp_params(ir);
2277 ir->fepvals->n_lambda = 0;
2280 /* WALL PARAMETERS */
2282 do_wall_params(ir, is->wall_atomtype, is->wall_density, opts);
2284 /* ORIENTATION RESTRAINT PARAMETERS */
2286 if (opts->bOrire && str_nelem(is->orirefitgrp, MAXPTR, NULL) != 1)
2288 warning_error(wi, "ERROR: Need one orientation restraint fit group\n");
2291 /* DEFORMATION PARAMETERS */
2293 clear_mat(ir->deform);
2294 for (i = 0; i < 6; i++)
2298 m = sscanf(is->deform, "%lf %lf %lf %lf %lf %lf",
2299 &(dumdub[0][0]), &(dumdub[0][1]), &(dumdub[0][2]),
2300 &(dumdub[0][3]), &(dumdub[0][4]), &(dumdub[0][5]));
2301 for (i = 0; i < 3; i++)
2303 ir->deform[i][i] = dumdub[0][i];
2305 ir->deform[YY][XX] = dumdub[0][3];
2306 ir->deform[ZZ][XX] = dumdub[0][4];
2307 ir->deform[ZZ][YY] = dumdub[0][5];
2308 if (ir->epc != epcNO)
2310 for (i = 0; i < 3; i++)
2312 for (j = 0; j <= i; j++)
2314 if (ir->deform[i][j] != 0 && ir->compress[i][j] != 0)
2316 warning_error(wi, "A box element has deform set and compressibility > 0");
2320 for (i = 0; i < 3; i++)
2322 for (j = 0; j < i; j++)
2324 if (ir->deform[i][j] != 0)
2326 for (m = j; m < DIM; m++)
2328 if (ir->compress[m][j] != 0)
2330 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.");
2331 warning(wi, warn_buf);
2343 static int search_QMstring(const char *s, int ng, const char *gn[])
2345 /* same as normal search_string, but this one searches QM strings */
2348 for (i = 0; (i < ng); i++)
2350 if (gmx_strcasecmp(s, gn[i]) == 0)
2356 gmx_fatal(FARGS, "this QM method or basisset (%s) is not implemented\n!", s);
2360 } /* search_QMstring */
2362 /* We would like gn to be const as well, but C doesn't allow this */
2363 int search_string(const char *s, int ng, char *gn[])
2367 for (i = 0; (i < ng); i++)
2369 if (gmx_strcasecmp(s, gn[i]) == 0)
2376 "Group %s referenced in the .mdp file was not found in the index file.\n"
2377 "Group names must match either [moleculetype] names or custom index group\n"
2378 "names, in which case you must supply an index file to the '-n' option\n"
2385 static gmx_bool do_numbering(int natoms, gmx_groups_t *groups, int ng, char *ptrs[],
2386 t_blocka *block, char *gnames[],
2387 int gtype, int restnm,
2388 int grptp, gmx_bool bVerbose,
2391 unsigned short *cbuf;
2392 t_grps *grps = &(groups->grps[gtype]);
2393 int i, j, gid, aj, ognr, ntot = 0;
2396 char warn_buf[STRLEN];
2400 fprintf(debug, "Starting numbering %d groups of type %d\n", ng, gtype);
2403 title = gtypes[gtype];
2406 /* Mark all id's as not set */
2407 for (i = 0; (i < natoms); i++)
2412 snew(grps->nm_ind, ng+1); /* +1 for possible rest group */
2413 for (i = 0; (i < ng); i++)
2415 /* Lookup the group name in the block structure */
2416 gid = search_string(ptrs[i], block->nr, gnames);
2417 if ((grptp != egrptpONE) || (i == 0))
2419 grps->nm_ind[grps->nr++] = gid;
2423 fprintf(debug, "Found gid %d for group %s\n", gid, ptrs[i]);
2426 /* Now go over the atoms in the group */
2427 for (j = block->index[gid]; (j < block->index[gid+1]); j++)
2432 /* Range checking */
2433 if ((aj < 0) || (aj >= natoms))
2435 gmx_fatal(FARGS, "Invalid atom number %d in indexfile", aj);
2437 /* Lookup up the old group number */
2441 gmx_fatal(FARGS, "Atom %d in multiple %s groups (%d and %d)",
2442 aj+1, title, ognr+1, i+1);
2446 /* Store the group number in buffer */
2447 if (grptp == egrptpONE)
2460 /* Now check whether we have done all atoms */
2464 if (grptp == egrptpALL)
2466 gmx_fatal(FARGS, "%d atoms are not part of any of the %s groups",
2467 natoms-ntot, title);
2469 else if (grptp == egrptpPART)
2471 sprintf(warn_buf, "%d atoms are not part of any of the %s groups",
2472 natoms-ntot, title);
2473 warning_note(wi, warn_buf);
2475 /* Assign all atoms currently unassigned to a rest group */
2476 for (j = 0; (j < natoms); j++)
2478 if (cbuf[j] == NOGID)
2484 if (grptp != egrptpPART)
2489 "Making dummy/rest group for %s containing %d elements\n",
2490 title, natoms-ntot);
2492 /* Add group name "rest" */
2493 grps->nm_ind[grps->nr] = restnm;
2495 /* Assign the rest name to all atoms not currently assigned to a group */
2496 for (j = 0; (j < natoms); j++)
2498 if (cbuf[j] == NOGID)
2507 if (grps->nr == 1 && (ntot == 0 || ntot == natoms))
2509 /* All atoms are part of one (or no) group, no index required */
2510 groups->ngrpnr[gtype] = 0;
2511 groups->grpnr[gtype] = NULL;
2515 groups->ngrpnr[gtype] = natoms;
2516 snew(groups->grpnr[gtype], natoms);
2517 for (j = 0; (j < natoms); j++)
2519 groups->grpnr[gtype][j] = cbuf[j];
2525 return (bRest && grptp == egrptpPART);
2528 static void calc_nrdf(gmx_mtop_t *mtop, t_inputrec *ir, char **gnames)
2531 gmx_groups_t *groups;
2533 int natoms, ai, aj, i, j, d, g, imin, jmin;
2535 int *nrdf2, *na_vcm, na_tot;
2536 double *nrdf_tc, *nrdf_vcm, nrdf_uc, n_sub = 0;
2537 gmx_mtop_atomloop_all_t aloop;
2539 int mb, mol, ftype, as;
2540 gmx_molblock_t *molb;
2541 gmx_moltype_t *molt;
2544 * First calc 3xnr-atoms for each group
2545 * then subtract half a degree of freedom for each constraint
2547 * Only atoms and nuclei contribute to the degrees of freedom...
2552 groups = &mtop->groups;
2553 natoms = mtop->natoms;
2555 /* Allocate one more for a possible rest group */
2556 /* We need to sum degrees of freedom into doubles,
2557 * since floats give too low nrdf's above 3 million atoms.
2559 snew(nrdf_tc, groups->grps[egcTC].nr+1);
2560 snew(nrdf_vcm, groups->grps[egcVCM].nr+1);
2561 snew(na_vcm, groups->grps[egcVCM].nr+1);
2563 for (i = 0; i < groups->grps[egcTC].nr; i++)
2567 for (i = 0; i < groups->grps[egcVCM].nr+1; i++)
2572 snew(nrdf2, natoms);
2573 aloop = gmx_mtop_atomloop_all_init(mtop);
2574 while (gmx_mtop_atomloop_all_next(aloop, &i, &atom))
2577 if (atom->ptype == eptAtom || atom->ptype == eptNucleus)
2579 g = ggrpnr(groups, egcFREEZE, i);
2580 /* Double count nrdf for particle i */
2581 for (d = 0; d < DIM; d++)
2583 if (opts->nFreeze[g][d] == 0)
2588 nrdf_tc [ggrpnr(groups, egcTC, i)] += 0.5*nrdf2[i];
2589 nrdf_vcm[ggrpnr(groups, egcVCM, i)] += 0.5*nrdf2[i];
2594 for (mb = 0; mb < mtop->nmolblock; mb++)
2596 molb = &mtop->molblock[mb];
2597 molt = &mtop->moltype[molb->type];
2598 atom = molt->atoms.atom;
2599 for (mol = 0; mol < molb->nmol; mol++)
2601 for (ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++)
2603 ia = molt->ilist[ftype].iatoms;
2604 for (i = 0; i < molt->ilist[ftype].nr; )
2606 /* Subtract degrees of freedom for the constraints,
2607 * if the particles still have degrees of freedom left.
2608 * If one of the particles is a vsite or a shell, then all
2609 * constraint motion will go there, but since they do not
2610 * contribute to the constraints the degrees of freedom do not
2615 if (((atom[ia[1]].ptype == eptNucleus) ||
2616 (atom[ia[1]].ptype == eptAtom)) &&
2617 ((atom[ia[2]].ptype == eptNucleus) ||
2618 (atom[ia[2]].ptype == eptAtom)))
2636 imin = min(imin, nrdf2[ai]);
2637 jmin = min(jmin, nrdf2[aj]);
2640 nrdf_tc [ggrpnr(groups, egcTC, ai)] -= 0.5*imin;
2641 nrdf_tc [ggrpnr(groups, egcTC, aj)] -= 0.5*jmin;
2642 nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5*imin;
2643 nrdf_vcm[ggrpnr(groups, egcVCM, aj)] -= 0.5*jmin;
2645 ia += interaction_function[ftype].nratoms+1;
2646 i += interaction_function[ftype].nratoms+1;
2649 ia = molt->ilist[F_SETTLE].iatoms;
2650 for (i = 0; i < molt->ilist[F_SETTLE].nr; )
2652 /* Subtract 1 dof from every atom in the SETTLE */
2653 for (j = 0; j < 3; j++)
2656 imin = min(2, nrdf2[ai]);
2658 nrdf_tc [ggrpnr(groups, egcTC, ai)] -= 0.5*imin;
2659 nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5*imin;
2664 as += molt->atoms.nr;
2668 if (ir->ePull == epullCONSTRAINT)
2670 /* Correct nrdf for the COM constraints.
2671 * We correct using the TC and VCM group of the first atom
2672 * in the reference and pull group. If atoms in one pull group
2673 * belong to different TC or VCM groups it is anyhow difficult
2674 * to determine the optimal nrdf assignment.
2678 for (i = 0; i < pull->ncoord; i++)
2682 for (j = 0; j < 2; j++)
2684 const t_pull_group *pgrp;
2686 pgrp = &pull->group[pull->coord[i].group[j]];
2690 /* Subtract 1/2 dof from each group */
2692 nrdf_tc [ggrpnr(groups, egcTC, ai)] -= 0.5*imin;
2693 nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5*imin;
2694 if (nrdf_tc[ggrpnr(groups, egcTC, ai)] < 0)
2696 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)]]);
2701 /* We need to subtract the whole DOF from group j=1 */
2708 if (ir->nstcomm != 0)
2710 /* Subtract 3 from the number of degrees of freedom in each vcm group
2711 * when com translation is removed and 6 when rotation is removed
2714 switch (ir->comm_mode)
2717 n_sub = ndof_com(ir);
2724 gmx_incons("Checking comm_mode");
2727 for (i = 0; i < groups->grps[egcTC].nr; i++)
2729 /* Count the number of atoms of TC group i for every VCM group */
2730 for (j = 0; j < groups->grps[egcVCM].nr+1; j++)
2735 for (ai = 0; ai < natoms; ai++)
2737 if (ggrpnr(groups, egcTC, ai) == i)
2739 na_vcm[ggrpnr(groups, egcVCM, ai)]++;
2743 /* Correct for VCM removal according to the fraction of each VCM
2744 * group present in this TC group.
2746 nrdf_uc = nrdf_tc[i];
2749 fprintf(debug, "T-group[%d] nrdf_uc = %g, n_sub = %g\n",
2753 for (j = 0; j < groups->grps[egcVCM].nr+1; j++)
2755 if (nrdf_vcm[j] > n_sub)
2757 nrdf_tc[i] += nrdf_uc*((double)na_vcm[j]/(double)na_tot)*
2758 (nrdf_vcm[j] - n_sub)/nrdf_vcm[j];
2762 fprintf(debug, " nrdf_vcm[%d] = %g, nrdf = %g\n",
2763 j, nrdf_vcm[j], nrdf_tc[i]);
2768 for (i = 0; (i < groups->grps[egcTC].nr); i++)
2770 opts->nrdf[i] = nrdf_tc[i];
2771 if (opts->nrdf[i] < 0)
2776 "Number of degrees of freedom in T-Coupling group %s is %.2f\n",
2777 gnames[groups->grps[egcTC].nm_ind[i]], opts->nrdf[i]);
2786 static void decode_cos(char *s, t_cosines *cosine)
2789 char format[STRLEN], f1[STRLEN];
2801 sscanf(t, "%d", &(cosine->n));
2808 snew(cosine->a, cosine->n);
2809 snew(cosine->phi, cosine->n);
2811 sprintf(format, "%%*d");
2812 for (i = 0; (i < cosine->n); i++)
2815 strcat(f1, "%lf%lf");
2816 if (sscanf(t, f1, &a, &phi) < 2)
2818 gmx_fatal(FARGS, "Invalid input for electric field shift: '%s'", t);
2821 cosine->phi[i] = phi;
2822 strcat(format, "%*lf%*lf");
2829 static gmx_bool do_egp_flag(t_inputrec *ir, gmx_groups_t *groups,
2830 const char *option, const char *val, int flag)
2832 /* The maximum number of energy group pairs would be MAXPTR*(MAXPTR+1)/2.
2833 * But since this is much larger than STRLEN, such a line can not be parsed.
2834 * The real maximum is the number of names that fit in a string: STRLEN/2.
2836 #define EGP_MAX (STRLEN/2)
2837 int nelem, i, j, k, nr;
2838 char *names[EGP_MAX];
2842 gnames = groups->grpname;
2844 nelem = str_nelem(val, EGP_MAX, names);
2847 gmx_fatal(FARGS, "The number of groups for %s is odd", option);
2849 nr = groups->grps[egcENER].nr;
2851 for (i = 0; i < nelem/2; i++)
2855 gmx_strcasecmp(names[2*i], *(gnames[groups->grps[egcENER].nm_ind[j]])))
2861 gmx_fatal(FARGS, "%s in %s is not an energy group\n",
2862 names[2*i], option);
2866 gmx_strcasecmp(names[2*i+1], *(gnames[groups->grps[egcENER].nm_ind[k]])))
2872 gmx_fatal(FARGS, "%s in %s is not an energy group\n",
2873 names[2*i+1], option);
2875 if ((j < nr) && (k < nr))
2877 ir->opts.egp_flags[nr*j+k] |= flag;
2878 ir->opts.egp_flags[nr*k+j] |= flag;
2886 void do_index(const char* mdparin, const char *ndx,
2889 t_inputrec *ir, rvec *v,
2893 gmx_groups_t *groups;
2897 char warnbuf[STRLEN], **gnames;
2898 int nr, ntcg, ntau_t, nref_t, nacc, nofg, nSA, nSA_points, nSA_time, nSA_temp;
2901 int nacg, nfreeze, nfrdim, nenergy, nvcm, nuser;
2902 char *ptr1[MAXPTR], *ptr2[MAXPTR], *ptr3[MAXPTR];
2903 int i, j, k, restnm;
2905 gmx_bool bExcl, bTable, bSetTCpar, bAnneal, bRest;
2906 int nQMmethod, nQMbasis, nQMcharge, nQMmult, nbSH, nCASorb, nCASelec,
2907 nSAon, nSAoff, nSAsteps, nQMg, nbOPT, nbTS;
2908 char warn_buf[STRLEN];
2912 fprintf(stderr, "processing index file...\n");
2918 snew(grps->index, 1);
2920 atoms_all = gmx_mtop_global_atoms(mtop);
2921 analyse(&atoms_all, grps, &gnames, FALSE, TRUE);
2922 free_t_atoms(&atoms_all, FALSE);
2926 grps = init_index(ndx, &gnames);
2929 groups = &mtop->groups;
2930 natoms = mtop->natoms;
2931 symtab = &mtop->symtab;
2933 snew(groups->grpname, grps->nr+1);
2935 for (i = 0; (i < grps->nr); i++)
2937 groups->grpname[i] = put_symtab(symtab, gnames[i]);
2939 groups->grpname[i] = put_symtab(symtab, "rest");
2941 srenew(gnames, grps->nr+1);
2942 gnames[restnm] = *(groups->grpname[i]);
2943 groups->ngrpname = grps->nr+1;
2945 set_warning_line(wi, mdparin, -1);
2947 ntau_t = str_nelem(is->tau_t, MAXPTR, ptr1);
2948 nref_t = str_nelem(is->ref_t, MAXPTR, ptr2);
2949 ntcg = str_nelem(is->tcgrps, MAXPTR, ptr3);
2950 if ((ntau_t != ntcg) || (nref_t != ntcg))
2952 gmx_fatal(FARGS, "Invalid T coupling input: %d groups, %d ref-t values and "
2953 "%d tau-t values", ntcg, nref_t, ntau_t);
2956 bSetTCpar = (ir->etc || EI_SD(ir->eI) || ir->eI == eiBD || EI_TPI(ir->eI));
2957 do_numbering(natoms, groups, ntcg, ptr3, grps, gnames, egcTC,
2958 restnm, bSetTCpar ? egrptpALL : egrptpALL_GENREST, bVerbose, wi);
2959 nr = groups->grps[egcTC].nr;
2961 snew(ir->opts.nrdf, nr);
2962 snew(ir->opts.tau_t, nr);
2963 snew(ir->opts.ref_t, nr);
2964 if (ir->eI == eiBD && ir->bd_fric == 0)
2966 fprintf(stderr, "bd-fric=0, so tau-t will be used as the inverse friction constant(s)\n");
2973 gmx_fatal(FARGS, "Not enough ref-t and tau-t values!");
2977 for (i = 0; (i < nr); i++)
2979 ir->opts.tau_t[i] = strtod(ptr1[i], NULL);
2980 if ((ir->eI == eiBD || ir->eI == eiSD2) && ir->opts.tau_t[i] <= 0)
2982 sprintf(warn_buf, "With integrator %s tau-t should be larger than 0", ei_names[ir->eI]);
2983 warning_error(wi, warn_buf);
2986 if (ir->etc != etcVRESCALE && ir->opts.tau_t[i] == 0)
2988 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.");
2991 if (ir->opts.tau_t[i] >= 0)
2993 tau_min = min(tau_min, ir->opts.tau_t[i]);
2996 if (ir->etc != etcNO && ir->nsttcouple == -1)
2998 ir->nsttcouple = ir_optimal_nsttcouple(ir);
3003 if ((ir->etc == etcNOSEHOOVER) && (ir->epc == epcBERENDSEN))
3005 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");
3007 if ((ir->epc == epcMTTK) && (ir->etc > etcNO))
3009 if (ir->nstpcouple != ir->nsttcouple)
3011 int mincouple = min(ir->nstpcouple, ir->nsttcouple);
3012 ir->nstpcouple = ir->nsttcouple = mincouple;
3013 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);
3014 warning_note(wi, warn_buf);
3018 /* velocity verlet with averaged kinetic energy KE = 0.5*(v(t+1/2) - v(t-1/2)) is implemented
3019 primarily for testing purposes, and does not work with temperature coupling other than 1 */
3021 if (ETC_ANDERSEN(ir->etc))
3023 if (ir->nsttcouple != 1)
3026 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");
3027 warning_note(wi, warn_buf);
3030 nstcmin = tcouple_min_integration_steps(ir->etc);
3033 if (tau_min/(ir->delta_t*ir->nsttcouple) < nstcmin)
3035 sprintf(warn_buf, "For proper integration of the %s thermostat, tau-t (%g) should be at least %d times larger than nsttcouple*dt (%g)",
3036 ETCOUPLTYPE(ir->etc),
3038 ir->nsttcouple*ir->delta_t);
3039 warning(wi, warn_buf);
3042 for (i = 0; (i < nr); i++)
3044 ir->opts.ref_t[i] = strtod(ptr2[i], NULL);
3045 if (ir->opts.ref_t[i] < 0)
3047 gmx_fatal(FARGS, "ref-t for group %d negative", i);
3050 /* set the lambda mc temperature to the md integrator temperature (which should be defined
3051 if we are in this conditional) if mc_temp is negative */
3052 if (ir->expandedvals->mc_temp < 0)
3054 ir->expandedvals->mc_temp = ir->opts.ref_t[0]; /*for now, set to the first reft */
3058 /* Simulated annealing for each group. There are nr groups */
3059 nSA = str_nelem(is->anneal, MAXPTR, ptr1);
3060 if (nSA == 1 && (ptr1[0][0] == 'n' || ptr1[0][0] == 'N'))
3064 if (nSA > 0 && nSA != nr)
3066 gmx_fatal(FARGS, "Not enough annealing values: %d (for %d groups)\n", nSA, nr);
3070 snew(ir->opts.annealing, nr);
3071 snew(ir->opts.anneal_npoints, nr);
3072 snew(ir->opts.anneal_time, nr);
3073 snew(ir->opts.anneal_temp, nr);
3074 for (i = 0; i < nr; i++)
3076 ir->opts.annealing[i] = eannNO;
3077 ir->opts.anneal_npoints[i] = 0;
3078 ir->opts.anneal_time[i] = NULL;
3079 ir->opts.anneal_temp[i] = NULL;
3084 for (i = 0; i < nr; i++)
3086 if (ptr1[i][0] == 'n' || ptr1[i][0] == 'N')
3088 ir->opts.annealing[i] = eannNO;
3090 else if (ptr1[i][0] == 's' || ptr1[i][0] == 'S')
3092 ir->opts.annealing[i] = eannSINGLE;
3095 else if (ptr1[i][0] == 'p' || ptr1[i][0] == 'P')
3097 ir->opts.annealing[i] = eannPERIODIC;
3103 /* Read the other fields too */
3104 nSA_points = str_nelem(is->anneal_npoints, MAXPTR, ptr1);
3105 if (nSA_points != nSA)
3107 gmx_fatal(FARGS, "Found %d annealing-npoints values for %d groups\n", nSA_points, nSA);
3109 for (k = 0, i = 0; i < nr; i++)
3111 ir->opts.anneal_npoints[i] = strtol(ptr1[i], NULL, 10);
3112 if (ir->opts.anneal_npoints[i] == 1)
3114 gmx_fatal(FARGS, "Please specify at least a start and an end point for annealing\n");
3116 snew(ir->opts.anneal_time[i], ir->opts.anneal_npoints[i]);
3117 snew(ir->opts.anneal_temp[i], ir->opts.anneal_npoints[i]);
3118 k += ir->opts.anneal_npoints[i];
3121 nSA_time = str_nelem(is->anneal_time, MAXPTR, ptr1);
3124 gmx_fatal(FARGS, "Found %d annealing-time values, wanter %d\n", nSA_time, k);
3126 nSA_temp = str_nelem(is->anneal_temp, MAXPTR, ptr2);
3129 gmx_fatal(FARGS, "Found %d annealing-temp values, wanted %d\n", nSA_temp, k);
3132 for (i = 0, k = 0; i < nr; i++)
3135 for (j = 0; j < ir->opts.anneal_npoints[i]; j++)
3137 ir->opts.anneal_time[i][j] = strtod(ptr1[k], NULL);
3138 ir->opts.anneal_temp[i][j] = strtod(ptr2[k], NULL);
3141 if (ir->opts.anneal_time[i][0] > (ir->init_t+GMX_REAL_EPS))
3143 gmx_fatal(FARGS, "First time point for annealing > init_t.\n");
3149 if (ir->opts.anneal_time[i][j] < ir->opts.anneal_time[i][j-1])
3151 gmx_fatal(FARGS, "Annealing timepoints out of order: t=%f comes after t=%f\n",
3152 ir->opts.anneal_time[i][j], ir->opts.anneal_time[i][j-1]);
3155 if (ir->opts.anneal_temp[i][j] < 0)
3157 gmx_fatal(FARGS, "Found negative temperature in annealing: %f\n", ir->opts.anneal_temp[i][j]);
3162 /* Print out some summary information, to make sure we got it right */
3163 for (i = 0, k = 0; i < nr; i++)
3165 if (ir->opts.annealing[i] != eannNO)
3167 j = groups->grps[egcTC].nm_ind[i];
3168 fprintf(stderr, "Simulated annealing for group %s: %s, %d timepoints\n",
3169 *(groups->grpname[j]), eann_names[ir->opts.annealing[i]],
3170 ir->opts.anneal_npoints[i]);
3171 fprintf(stderr, "Time (ps) Temperature (K)\n");
3172 /* All terms except the last one */
3173 for (j = 0; j < (ir->opts.anneal_npoints[i]-1); j++)
3175 fprintf(stderr, "%9.1f %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3178 /* Finally the last one */
3179 j = ir->opts.anneal_npoints[i]-1;
3180 if (ir->opts.annealing[i] == eannSINGLE)
3182 fprintf(stderr, "%9.1f- %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3186 fprintf(stderr, "%9.1f %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3187 if (fabs(ir->opts.anneal_temp[i][j]-ir->opts.anneal_temp[i][0]) > GMX_REAL_EPS)
3189 warning_note(wi, "There is a temperature jump when your annealing loops back.\n");
3198 if (ir->ePull != epullNO)
3200 make_pull_groups(ir->pull, is->pull_grp, grps, gnames);
3202 make_pull_coords(ir->pull);
3207 make_rotation_groups(ir->rot, is->rot_grp, grps, gnames);
3210 nacc = str_nelem(is->acc, MAXPTR, ptr1);
3211 nacg = str_nelem(is->accgrps, MAXPTR, ptr2);
3212 if (nacg*DIM != nacc)
3214 gmx_fatal(FARGS, "Invalid Acceleration input: %d groups and %d acc. values",
3217 do_numbering(natoms, groups, nacg, ptr2, grps, gnames, egcACC,
3218 restnm, egrptpALL_GENREST, bVerbose, wi);
3219 nr = groups->grps[egcACC].nr;
3220 snew(ir->opts.acc, nr);
3221 ir->opts.ngacc = nr;
3223 for (i = k = 0; (i < nacg); i++)
3225 for (j = 0; (j < DIM); j++, k++)
3227 ir->opts.acc[i][j] = strtod(ptr1[k], NULL);
3230 for (; (i < nr); i++)
3232 for (j = 0; (j < DIM); j++)
3234 ir->opts.acc[i][j] = 0;
3238 nfrdim = str_nelem(is->frdim, MAXPTR, ptr1);
3239 nfreeze = str_nelem(is->freeze, MAXPTR, ptr2);
3240 if (nfrdim != DIM*nfreeze)
3242 gmx_fatal(FARGS, "Invalid Freezing input: %d groups and %d freeze values",
3245 do_numbering(natoms, groups, nfreeze, ptr2, grps, gnames, egcFREEZE,
3246 restnm, egrptpALL_GENREST, bVerbose, wi);
3247 nr = groups->grps[egcFREEZE].nr;
3248 ir->opts.ngfrz = nr;
3249 snew(ir->opts.nFreeze, nr);
3250 for (i = k = 0; (i < nfreeze); i++)
3252 for (j = 0; (j < DIM); j++, k++)
3254 ir->opts.nFreeze[i][j] = (gmx_strncasecmp(ptr1[k], "Y", 1) == 0);
3255 if (!ir->opts.nFreeze[i][j])
3257 if (gmx_strncasecmp(ptr1[k], "N", 1) != 0)
3259 sprintf(warnbuf, "Please use Y(ES) or N(O) for freezedim only "
3260 "(not %s)", ptr1[k]);
3261 warning(wi, warn_buf);
3266 for (; (i < nr); i++)
3268 for (j = 0; (j < DIM); j++)
3270 ir->opts.nFreeze[i][j] = 0;
3274 nenergy = str_nelem(is->energy, MAXPTR, ptr1);
3275 do_numbering(natoms, groups, nenergy, ptr1, grps, gnames, egcENER,
3276 restnm, egrptpALL_GENREST, bVerbose, wi);
3277 add_wall_energrps(groups, ir->nwall, symtab);
3278 ir->opts.ngener = groups->grps[egcENER].nr;
3279 nvcm = str_nelem(is->vcm, MAXPTR, ptr1);
3281 do_numbering(natoms, groups, nvcm, ptr1, grps, gnames, egcVCM,
3282 restnm, nvcm == 0 ? egrptpALL_GENREST : egrptpPART, bVerbose, wi);
3285 warning(wi, "Some atoms are not part of any center of mass motion removal group.\n"
3286 "This may lead to artifacts.\n"
3287 "In most cases one should use one group for the whole system.");
3290 /* Now we have filled the freeze struct, so we can calculate NRDF */
3291 calc_nrdf(mtop, ir, gnames);
3297 /* Must check per group! */
3298 for (i = 0; (i < ir->opts.ngtc); i++)
3300 ntot += ir->opts.nrdf[i];
3302 if (ntot != (DIM*natoms))
3304 fac = sqrt(ntot/(DIM*natoms));
3307 fprintf(stderr, "Scaling velocities by a factor of %.3f to account for constraints\n"
3308 "and removal of center of mass motion\n", fac);
3310 for (i = 0; (i < natoms); i++)
3312 svmul(fac, v[i], v[i]);
3317 nuser = str_nelem(is->user1, MAXPTR, ptr1);
3318 do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcUser1,
3319 restnm, egrptpALL_GENREST, bVerbose, wi);
3320 nuser = str_nelem(is->user2, MAXPTR, ptr1);
3321 do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcUser2,
3322 restnm, egrptpALL_GENREST, bVerbose, wi);
3323 nuser = str_nelem(is->xtc_grps, MAXPTR, ptr1);
3324 do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcXTC,
3325 restnm, egrptpONE, bVerbose, wi);
3326 nofg = str_nelem(is->orirefitgrp, MAXPTR, ptr1);
3327 do_numbering(natoms, groups, nofg, ptr1, grps, gnames, egcORFIT,
3328 restnm, egrptpALL_GENREST, bVerbose, wi);
3330 /* QMMM input processing */
3331 nQMg = str_nelem(is->QMMM, MAXPTR, ptr1);
3332 nQMmethod = str_nelem(is->QMmethod, MAXPTR, ptr2);
3333 nQMbasis = str_nelem(is->QMbasis, MAXPTR, ptr3);
3334 if ((nQMmethod != nQMg) || (nQMbasis != nQMg))
3336 gmx_fatal(FARGS, "Invalid QMMM input: %d groups %d basissets"
3337 " and %d methods\n", nQMg, nQMbasis, nQMmethod);
3339 /* group rest, if any, is always MM! */
3340 do_numbering(natoms, groups, nQMg, ptr1, grps, gnames, egcQMMM,
3341 restnm, egrptpALL_GENREST, bVerbose, wi);
3342 nr = nQMg; /*atoms->grps[egcQMMM].nr;*/
3343 ir->opts.ngQM = nQMg;
3344 snew(ir->opts.QMmethod, nr);
3345 snew(ir->opts.QMbasis, nr);
3346 for (i = 0; i < nr; i++)
3348 /* input consists of strings: RHF CASSCF PM3 .. These need to be
3349 * converted to the corresponding enum in names.c
3351 ir->opts.QMmethod[i] = search_QMstring(ptr2[i], eQMmethodNR,
3353 ir->opts.QMbasis[i] = search_QMstring(ptr3[i], eQMbasisNR,
3357 nQMmult = str_nelem(is->QMmult, MAXPTR, ptr1);
3358 nQMcharge = str_nelem(is->QMcharge, MAXPTR, ptr2);
3359 nbSH = str_nelem(is->bSH, MAXPTR, ptr3);
3360 snew(ir->opts.QMmult, nr);
3361 snew(ir->opts.QMcharge, nr);
3362 snew(ir->opts.bSH, nr);
3364 for (i = 0; i < nr; i++)
3366 ir->opts.QMmult[i] = strtol(ptr1[i], NULL, 10);
3367 ir->opts.QMcharge[i] = strtol(ptr2[i], NULL, 10);
3368 ir->opts.bSH[i] = (gmx_strncasecmp(ptr3[i], "Y", 1) == 0);
3371 nCASelec = str_nelem(is->CASelectrons, MAXPTR, ptr1);
3372 nCASorb = str_nelem(is->CASorbitals, MAXPTR, ptr2);
3373 snew(ir->opts.CASelectrons, nr);
3374 snew(ir->opts.CASorbitals, nr);
3375 for (i = 0; i < nr; i++)
3377 ir->opts.CASelectrons[i] = strtol(ptr1[i], NULL, 10);
3378 ir->opts.CASorbitals[i] = strtol(ptr2[i], NULL, 10);
3380 /* special optimization options */
3382 nbOPT = str_nelem(is->bOPT, MAXPTR, ptr1);
3383 nbTS = str_nelem(is->bTS, MAXPTR, ptr2);
3384 snew(ir->opts.bOPT, nr);
3385 snew(ir->opts.bTS, nr);
3386 for (i = 0; i < nr; i++)
3388 ir->opts.bOPT[i] = (gmx_strncasecmp(ptr1[i], "Y", 1) == 0);
3389 ir->opts.bTS[i] = (gmx_strncasecmp(ptr2[i], "Y", 1) == 0);
3391 nSAon = str_nelem(is->SAon, MAXPTR, ptr1);
3392 nSAoff = str_nelem(is->SAoff, MAXPTR, ptr2);
3393 nSAsteps = str_nelem(is->SAsteps, MAXPTR, ptr3);
3394 snew(ir->opts.SAon, nr);
3395 snew(ir->opts.SAoff, nr);
3396 snew(ir->opts.SAsteps, nr);
3398 for (i = 0; i < nr; i++)
3400 ir->opts.SAon[i] = strtod(ptr1[i], NULL);
3401 ir->opts.SAoff[i] = strtod(ptr2[i], NULL);
3402 ir->opts.SAsteps[i] = strtol(ptr3[i], NULL, 10);
3404 /* end of QMMM input */
3408 for (i = 0; (i < egcNR); i++)
3410 fprintf(stderr, "%-16s has %d element(s):", gtypes[i], groups->grps[i].nr);
3411 for (j = 0; (j < groups->grps[i].nr); j++)
3413 fprintf(stderr, " %s", *(groups->grpname[groups->grps[i].nm_ind[j]]));
3415 fprintf(stderr, "\n");
3419 nr = groups->grps[egcENER].nr;
3420 snew(ir->opts.egp_flags, nr*nr);
3422 bExcl = do_egp_flag(ir, groups, "energygrp-excl", is->egpexcl, EGP_EXCL);
3423 if (bExcl && ir->cutoff_scheme == ecutsVERLET)
3425 warning_error(wi, "Energy group exclusions are not (yet) implemented for the Verlet scheme");
3427 if (bExcl && EEL_FULL(ir->coulombtype))
3429 warning(wi, "Can not exclude the lattice Coulomb energy between energy groups");
3432 bTable = do_egp_flag(ir, groups, "energygrp-table", is->egptable, EGP_TABLE);
3433 if (bTable && !(ir->vdwtype == evdwUSER) &&
3434 !(ir->coulombtype == eelUSER) && !(ir->coulombtype == eelPMEUSER) &&
3435 !(ir->coulombtype == eelPMEUSERSWITCH))
3437 gmx_fatal(FARGS, "Can only have energy group pair tables in combination with user tables for VdW and/or Coulomb");
3440 decode_cos(is->efield_x, &(ir->ex[XX]));
3441 decode_cos(is->efield_xt, &(ir->et[XX]));
3442 decode_cos(is->efield_y, &(ir->ex[YY]));
3443 decode_cos(is->efield_yt, &(ir->et[YY]));
3444 decode_cos(is->efield_z, &(ir->ex[ZZ]));
3445 decode_cos(is->efield_zt, &(ir->et[ZZ]));
3449 do_adress_index(ir->adress, groups, gnames, &(ir->opts), wi);
3452 for (i = 0; (i < grps->nr); i++)
3464 static void check_disre(gmx_mtop_t *mtop)
3466 gmx_ffparams_t *ffparams;
3467 t_functype *functype;
3469 int i, ndouble, ftype;
3470 int label, old_label;
3472 if (gmx_mtop_ftype_count(mtop, F_DISRES) > 0)
3474 ffparams = &mtop->ffparams;
3475 functype = ffparams->functype;
3476 ip = ffparams->iparams;
3479 for (i = 0; i < ffparams->ntypes; i++)
3481 ftype = functype[i];
3482 if (ftype == F_DISRES)
3484 label = ip[i].disres.label;
3485 if (label == old_label)
3487 fprintf(stderr, "Distance restraint index %d occurs twice\n", label);
3495 gmx_fatal(FARGS, "Found %d double distance restraint indices,\n"
3496 "probably the parameters for multiple pairs in one restraint "
3497 "are not identical\n", ndouble);
3502 static gmx_bool absolute_reference(t_inputrec *ir, gmx_mtop_t *sys,
3503 gmx_bool posres_only,
3507 gmx_mtop_ilistloop_t iloop;
3517 for (d = 0; d < DIM; d++)
3519 AbsRef[d] = (d < ndof_com(ir) ? 0 : 1);
3521 /* Check for freeze groups */
3522 for (g = 0; g < ir->opts.ngfrz; g++)
3524 for (d = 0; d < DIM; d++)
3526 if (ir->opts.nFreeze[g][d] != 0)
3534 /* Check for position restraints */
3535 iloop = gmx_mtop_ilistloop_init(sys);
3536 while (gmx_mtop_ilistloop_next(iloop, &ilist, &nmol))
3539 (AbsRef[XX] == 0 || AbsRef[YY] == 0 || AbsRef[ZZ] == 0))
3541 for (i = 0; i < ilist[F_POSRES].nr; i += 2)
3543 pr = &sys->ffparams.iparams[ilist[F_POSRES].iatoms[i]];
3544 for (d = 0; d < DIM; d++)
3546 if (pr->posres.fcA[d] != 0)
3552 for (i = 0; i < ilist[F_FBPOSRES].nr; i += 2)
3554 /* Check for flat-bottom posres */
3555 pr = &sys->ffparams.iparams[ilist[F_FBPOSRES].iatoms[i]];
3556 if (pr->fbposres.k != 0)
3558 switch (pr->fbposres.geom)
3560 case efbposresSPHERE:
3561 AbsRef[XX] = AbsRef[YY] = AbsRef[ZZ] = 1;
3563 case efbposresCYLINDER:
3564 AbsRef[XX] = AbsRef[YY] = 1;
3566 case efbposresX: /* d=XX */
3567 case efbposresY: /* d=YY */
3568 case efbposresZ: /* d=ZZ */
3569 d = pr->fbposres.geom - efbposresX;
3573 gmx_fatal(FARGS, " Invalid geometry for flat-bottom position restraint.\n"
3574 "Expected nr between 1 and %d. Found %d\n", efbposresNR-1,
3582 return (AbsRef[XX] != 0 && AbsRef[YY] != 0 && AbsRef[ZZ] != 0);
3586 check_combination_rule_differences(const gmx_mtop_t *mtop, int state,
3587 gmx_bool *bC6ParametersWorkWithGeometricRules,
3588 gmx_bool *bC6ParametersWorkWithLBRules,
3589 gmx_bool *bLBRulesPossible)
3591 int ntypes, tpi, tpj, thisLBdiff, thisgeomdiff;
3594 double geometricdiff, LBdiff;
3595 double c6i, c6j, c12i, c12j;
3596 double c6, c6_geometric, c6_LB;
3597 double sigmai, sigmaj, epsi, epsj;
3598 gmx_bool bCanDoLBRules, bCanDoGeometricRules;
3601 /* A tolerance of 1e-5 seems reasonable for (possibly hand-typed)
3602 * force-field floating point parameters.
3605 ptr = getenv("GMX_LJCOMB_TOL");
3610 sscanf(ptr, "%lf", &dbl);
3614 *bC6ParametersWorkWithLBRules = TRUE;
3615 *bC6ParametersWorkWithGeometricRules = TRUE;
3616 bCanDoLBRules = TRUE;
3617 bCanDoGeometricRules = TRUE;
3618 ntypes = mtop->ffparams.atnr;
3619 snew(typecount, ntypes);
3620 gmx_mtop_count_atomtypes(mtop, state, typecount);
3621 geometricdiff = LBdiff = 0.0;
3622 *bLBRulesPossible = TRUE;
3623 for (tpi = 0; tpi < ntypes; ++tpi)
3625 c6i = mtop->ffparams.iparams[(ntypes + 1) * tpi].lj.c6;
3626 c12i = mtop->ffparams.iparams[(ntypes + 1) * tpi].lj.c12;
3627 for (tpj = tpi; tpj < ntypes; ++tpj)
3629 c6j = mtop->ffparams.iparams[(ntypes + 1) * tpj].lj.c6;
3630 c12j = mtop->ffparams.iparams[(ntypes + 1) * tpj].lj.c12;
3631 c6 = mtop->ffparams.iparams[ntypes * tpi + tpj].lj.c6;
3632 c6_geometric = sqrt(c6i * c6j);
3633 if (!gmx_numzero(c6_geometric))
3635 if (!gmx_numzero(c12i) && !gmx_numzero(c12j))
3637 sigmai = pow(c12i / c6i, 1.0/6.0);
3638 sigmaj = pow(c12j / c6j, 1.0/6.0);
3639 epsi = c6i * c6i /(4.0 * c12i);
3640 epsj = c6j * c6j /(4.0 * c12j);
3641 c6_LB = 4.0 * pow(epsi * epsj, 1.0/2.0) * pow(0.5 * (sigmai + sigmaj), 6);
3645 *bLBRulesPossible = FALSE;
3646 c6_LB = c6_geometric;
3648 bCanDoLBRules = gmx_within_tol(c6_LB, c6, tol);
3651 if (FALSE == bCanDoLBRules)
3653 *bC6ParametersWorkWithLBRules = FALSE;
3656 bCanDoGeometricRules = gmx_within_tol(c6_geometric, c6, tol);
3658 if (FALSE == bCanDoGeometricRules)
3660 *bC6ParametersWorkWithGeometricRules = FALSE;
3668 check_combination_rules(const t_inputrec *ir, const gmx_mtop_t *mtop,
3672 gmx_bool bLBRulesPossible, bC6ParametersWorkWithGeometricRules, bC6ParametersWorkWithLBRules;
3674 check_combination_rule_differences(mtop, 0,
3675 &bC6ParametersWorkWithGeometricRules,
3676 &bC6ParametersWorkWithLBRules,
3678 if (ir->ljpme_combination_rule == eljpmeLB)
3680 if (FALSE == bC6ParametersWorkWithLBRules || FALSE == bLBRulesPossible)
3682 warning(wi, "You are using arithmetic-geometric combination rules "
3683 "in LJ-PME, but your non-bonded C6 parameters do not "
3684 "follow these rules.");
3689 if (FALSE == bC6ParametersWorkWithGeometricRules)
3691 if (ir->eDispCorr != edispcNO)
3693 warning_note(wi, "You are using geometric combination rules in "
3694 "LJ-PME, but your non-bonded C6 parameters do "
3695 "not follow these rules. "
3696 "If your force field uses Lorentz-Berthelot combination rules this "
3697 "will introduce small errors in the forces and energies in "
3698 "your simulations. Dispersion correction will correct total "
3699 "energy and/or pressure, but not forces or surface tensions. "
3700 "Please check the LJ-PME section in the manual "
3701 "before proceeding further.");
3705 warning_note(wi, "You are using geometric combination rules in "
3706 "LJ-PME, but your non-bonded C6 parameters do "
3707 "not follow these rules. "
3708 "If your force field uses Lorentz-Berthelot combination rules this "
3709 "will introduce small errors in the forces and energies in "
3710 "your simulations. Consider using dispersion correction "
3711 "for the total energy and pressure. "
3712 "Please check the LJ-PME section in the manual "
3713 "before proceeding further.");
3719 void triple_check(const char *mdparin, t_inputrec *ir, gmx_mtop_t *sys,
3723 int i, m, c, nmol, npct;
3724 gmx_bool bCharge, bAcc;
3725 real gdt_max, *mgrp, mt;
3727 gmx_mtop_atomloop_block_t aloopb;
3728 gmx_mtop_atomloop_all_t aloop;
3731 char warn_buf[STRLEN];
3733 set_warning_line(wi, mdparin, -1);
3735 if (EI_DYNAMICS(ir->eI) && !EI_SD(ir->eI) && ir->eI != eiBD &&
3736 ir->comm_mode == ecmNO &&
3737 !(absolute_reference(ir, sys, FALSE, AbsRef) || ir->nsteps <= 10))
3739 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");
3742 /* Check for pressure coupling with absolute position restraints */
3743 if (ir->epc != epcNO && ir->refcoord_scaling == erscNO)
3745 absolute_reference(ir, sys, TRUE, AbsRef);
3747 for (m = 0; m < DIM; m++)
3749 if (AbsRef[m] && norm2(ir->compress[m]) > 0)
3751 warning(wi, "You are using pressure coupling with absolute position restraints, this will give artifacts. Use the refcoord_scaling option.");
3759 aloopb = gmx_mtop_atomloop_block_init(sys);
3760 while (gmx_mtop_atomloop_block_next(aloopb, &atom, &nmol))
3762 if (atom->q != 0 || atom->qB != 0)
3770 if (EEL_FULL(ir->coulombtype))
3773 "You are using full electrostatics treatment %s for a system without charges.\n"
3774 "This costs a lot of performance for just processing zeros, consider using %s instead.\n",
3775 EELTYPE(ir->coulombtype), EELTYPE(eelCUT));
3776 warning(wi, err_buf);
3781 if (ir->coulombtype == eelCUT && ir->rcoulomb > 0 && !ir->implicit_solvent)
3784 "You are using a plain Coulomb cut-off, which might produce artifacts.\n"
3785 "You might want to consider using %s electrostatics.\n",
3787 warning_note(wi, err_buf);
3791 /* Check if combination rules used in LJ-PME are the same as in the force field */
3792 if (EVDW_PME(ir->vdwtype))
3794 check_combination_rules(ir, sys, wi);
3797 /* Generalized reaction field */
3798 if (ir->opts.ngtc == 0)
3800 sprintf(err_buf, "No temperature coupling while using coulombtype %s",
3802 CHECK(ir->coulombtype == eelGRF);
3806 sprintf(err_buf, "When using coulombtype = %s"
3807 " ref-t for temperature coupling should be > 0",
3809 CHECK((ir->coulombtype == eelGRF) && (ir->opts.ref_t[0] <= 0));
3812 if (ir->eI == eiSD1 &&
3813 (gmx_mtop_ftype_count(sys, F_CONSTR) > 0 ||
3814 gmx_mtop_ftype_count(sys, F_SETTLE) > 0))
3816 sprintf(warn_buf, "With constraints integrator %s is less accurate, consider using %s instead", ei_names[ir->eI], ei_names[eiSD2]);
3817 warning_note(wi, warn_buf);
3821 for (i = 0; (i < sys->groups.grps[egcACC].nr); i++)
3823 for (m = 0; (m < DIM); m++)
3825 if (fabs(ir->opts.acc[i][m]) > 1e-6)
3834 snew(mgrp, sys->groups.grps[egcACC].nr);
3835 aloop = gmx_mtop_atomloop_all_init(sys);
3836 while (gmx_mtop_atomloop_all_next(aloop, &i, &atom))
3838 mgrp[ggrpnr(&sys->groups, egcACC, i)] += atom->m;
3841 for (i = 0; (i < sys->groups.grps[egcACC].nr); i++)
3843 for (m = 0; (m < DIM); m++)
3845 acc[m] += ir->opts.acc[i][m]*mgrp[i];
3849 for (m = 0; (m < DIM); m++)
3851 if (fabs(acc[m]) > 1e-6)
3853 const char *dim[DIM] = { "X", "Y", "Z" };
3855 "Net Acceleration in %s direction, will %s be corrected\n",
3856 dim[m], ir->nstcomm != 0 ? "" : "not");
3857 if (ir->nstcomm != 0 && m < ndof_com(ir))
3860 for (i = 0; (i < sys->groups.grps[egcACC].nr); i++)
3862 ir->opts.acc[i][m] -= acc[m];
3870 if (ir->efep != efepNO && ir->fepvals->sc_alpha != 0 &&
3871 !gmx_within_tol(sys->ffparams.reppow, 12.0, 10*GMX_DOUBLE_EPS))
3873 gmx_fatal(FARGS, "Soft-core interactions are only supported with VdW repulsion power 12");
3876 if (ir->ePull != epullNO)
3878 gmx_bool bPullAbsoluteRef;
3880 bPullAbsoluteRef = FALSE;
3881 for (i = 0; i < ir->pull->ncoord; i++)
3883 bPullAbsoluteRef = bPullAbsoluteRef ||
3884 ir->pull->coord[i].group[0] == 0 ||
3885 ir->pull->coord[i].group[1] == 0;
3887 if (bPullAbsoluteRef)
3889 absolute_reference(ir, sys, FALSE, AbsRef);
3890 for (m = 0; m < DIM; m++)
3892 if (ir->pull->dim[m] && !AbsRef[m])
3894 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.");
3900 if (ir->pull->eGeom == epullgDIRPBC)
3902 for (i = 0; i < 3; i++)
3904 for (m = 0; m <= i; m++)
3906 if ((ir->epc != epcNO && ir->compress[i][m] != 0) ||
3907 ir->deform[i][m] != 0)
3909 for (c = 0; c < ir->pull->ncoord; c++)
3911 if (ir->pull->coord[c].vec[m] != 0)
3913 gmx_fatal(FARGS, "Can not have dynamic box while using pull geometry '%s' (dim %c)", EPULLGEOM(ir->pull->eGeom), 'x'+m);
3925 void double_check(t_inputrec *ir, matrix box, gmx_bool bConstr, warninp_t wi)
3929 char warn_buf[STRLEN];
3932 ptr = check_box(ir->ePBC, box);
3935 warning_error(wi, ptr);
3938 if (bConstr && ir->eConstrAlg == econtSHAKE)
3940 if (ir->shake_tol <= 0.0)
3942 sprintf(warn_buf, "ERROR: shake-tol must be > 0 instead of %g\n",
3944 warning_error(wi, warn_buf);
3947 if (IR_TWINRANGE(*ir) && ir->nstlist > 1)
3949 sprintf(warn_buf, "With twin-range cut-off's and SHAKE the virial and the pressure are incorrect.");
3950 if (ir->epc == epcNO)
3952 warning(wi, warn_buf);
3956 warning_error(wi, warn_buf);
3961 if ( (ir->eConstrAlg == econtLINCS) && bConstr)
3963 /* If we have Lincs constraints: */
3964 if (ir->eI == eiMD && ir->etc == etcNO &&
3965 ir->eConstrAlg == econtLINCS && ir->nLincsIter == 1)
3967 sprintf(warn_buf, "For energy conservation with LINCS, lincs_iter should be 2 or larger.\n");
3968 warning_note(wi, warn_buf);
3971 if ((ir->eI == eiCG || ir->eI == eiLBFGS) && (ir->nProjOrder < 8))
3973 sprintf(warn_buf, "For accurate %s with LINCS constraints, lincs-order should be 8 or more.", ei_names[ir->eI]);
3974 warning_note(wi, warn_buf);
3976 if (ir->epc == epcMTTK)
3978 warning_error(wi, "MTTK not compatible with lincs -- use shake instead.");
3982 if (ir->LincsWarnAngle > 90.0)
3984 sprintf(warn_buf, "lincs-warnangle can not be larger than 90 degrees, setting it to 90.\n");
3985 warning(wi, warn_buf);
3986 ir->LincsWarnAngle = 90.0;
3989 if (ir->ePBC != epbcNONE)
3991 if (ir->nstlist == 0)
3993 warning(wi, "With nstlist=0 atoms are only put into the box at step 0, therefore drifting atoms might cause the simulation to crash.");
3995 bTWIN = (ir->rlistlong > ir->rlist);
3996 if (ir->ns_type == ensGRID)
3998 if (sqr(ir->rlistlong) >= max_cutoff2(ir->ePBC, box))
4000 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",
4001 bTWIN ? (ir->rcoulomb == ir->rlistlong ? "rcoulomb" : "rvdw") : "rlist");
4002 warning_error(wi, warn_buf);
4007 min_size = min(box[XX][XX], min(box[YY][YY], box[ZZ][ZZ]));
4008 if (2*ir->rlistlong >= min_size)
4010 sprintf(warn_buf, "ERROR: One of the box lengths is smaller than twice the cut-off length. Increase the box size or decrease rlist.");
4011 warning_error(wi, warn_buf);
4014 fprintf(stderr, "Grid search might allow larger cut-off's than simple search with triclinic boxes.");
4021 void check_chargegroup_radii(const gmx_mtop_t *mtop, const t_inputrec *ir,
4025 real rvdw1, rvdw2, rcoul1, rcoul2;
4026 char warn_buf[STRLEN];
4028 calc_chargegroup_radii(mtop, x, &rvdw1, &rvdw2, &rcoul1, &rcoul2);
4032 printf("Largest charge group radii for Van der Waals: %5.3f, %5.3f nm\n",
4037 printf("Largest charge group radii for Coulomb: %5.3f, %5.3f nm\n",
4043 if (rvdw1 + rvdw2 > ir->rlist ||
4044 rcoul1 + rcoul2 > ir->rlist)
4047 "The sum of the two largest charge group radii (%f) "
4048 "is larger than rlist (%f)\n",
4049 max(rvdw1+rvdw2, rcoul1+rcoul2), ir->rlist);
4050 warning(wi, warn_buf);
4054 /* Here we do not use the zero at cut-off macro,
4055 * since user defined interactions might purposely
4056 * not be zero at the cut-off.
4058 if ((EVDW_IS_ZERO_AT_CUTOFF(ir->vdwtype) ||
4059 ir->vdw_modifier != eintmodNONE) &&
4060 rvdw1 + rvdw2 > ir->rlistlong - ir->rvdw)
4062 sprintf(warn_buf, "The sum of the two largest charge group "
4063 "radii (%f) is larger than %s (%f) - rvdw (%f).\n"
4064 "With exact cut-offs, better performance can be "
4065 "obtained with cutoff-scheme = %s, because it "
4066 "does not use charge groups at all.",
4068 ir->rlistlong > ir->rlist ? "rlistlong" : "rlist",
4069 ir->rlistlong, ir->rvdw,
4070 ecutscheme_names[ecutsVERLET]);
4073 warning(wi, warn_buf);
4077 warning_note(wi, warn_buf);
4080 if ((EEL_IS_ZERO_AT_CUTOFF(ir->coulombtype) ||
4081 ir->coulomb_modifier != eintmodNONE) &&
4082 rcoul1 + rcoul2 > ir->rlistlong - ir->rcoulomb)
4084 sprintf(warn_buf, "The sum of the two largest charge group radii (%f) is larger than %s (%f) - rcoulomb (%f).\n"
4085 "With exact cut-offs, better performance can be obtained with cutoff-scheme = %s, because it does not use charge groups at all.",
4087 ir->rlistlong > ir->rlist ? "rlistlong" : "rlist",
4088 ir->rlistlong, ir->rcoulomb,
4089 ecutscheme_names[ecutsVERLET]);
4092 warning(wi, warn_buf);
4096 warning_note(wi, warn_buf);