2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
49 #include "gmx_fatal.h"
62 #include "mtop_util.h"
63 #include "chargegroup.h"
68 #define MAXLAMBDAS 1024
70 /* Resource parameters
71 * Do not change any of these until you read the instruction
72 * in readinp.h. Some cpp's do not take spaces after the backslash
73 * (like the c-shell), which will give you a very weird compiler
77 static char tcgrps[STRLEN], tau_t[STRLEN], ref_t[STRLEN],
78 acc[STRLEN], accgrps[STRLEN], freeze[STRLEN], frdim[STRLEN],
79 energy[STRLEN], user1[STRLEN], user2[STRLEN], vcm[STRLEN], xtc_grps[STRLEN],
80 couple_moltype[STRLEN], orirefitgrp[STRLEN], egptable[STRLEN], egpexcl[STRLEN],
81 wall_atomtype[STRLEN], wall_density[STRLEN], deform[STRLEN], QMMM[STRLEN];
82 static char fep_lambda[efptNR][STRLEN];
83 static char lambda_weights[STRLEN];
84 static char **pull_grp;
85 static char **rot_grp;
86 static char anneal[STRLEN], anneal_npoints[STRLEN],
87 anneal_time[STRLEN], anneal_temp[STRLEN];
88 static char QMmethod[STRLEN], QMbasis[STRLEN], QMcharge[STRLEN], QMmult[STRLEN],
89 bSH[STRLEN], CASorbitals[STRLEN], CASelectrons[STRLEN], SAon[STRLEN],
90 SAoff[STRLEN], SAsteps[STRLEN], bTS[STRLEN], bOPT[STRLEN];
91 static char efield_x[STRLEN], efield_xt[STRLEN], efield_y[STRLEN],
92 efield_yt[STRLEN], efield_z[STRLEN], efield_zt[STRLEN];
95 egrptpALL, /* All particles have to be a member of a group. */
96 egrptpALL_GENREST, /* A rest group with name is generated for particles *
97 * that are not part of any group. */
98 egrptpPART, /* As egrptpALL_GENREST, but no name is generated *
99 * for the rest group. */
100 egrptpONE /* Merge all selected groups into one group, *
101 * make a rest group for the remaining particles. */
104 static const char *constraints[eshNR+1] = {
105 "none", "h-bonds", "all-bonds", "h-angles", "all-angles", NULL
108 static const char *couple_lam[ecouplamNR+1] = {
109 "vdw-q", "vdw", "q", "none", NULL
112 void init_ir(t_inputrec *ir, t_gromppopts *opts)
114 snew(opts->include, STRLEN);
115 snew(opts->define, STRLEN);
116 snew(ir->fepvals, 1);
117 snew(ir->expandedvals, 1);
118 snew(ir->simtempvals, 1);
121 static void GetSimTemps(int ntemps, t_simtemp *simtemp, double *temperature_lambdas)
126 for (i = 0; i < ntemps; i++)
128 /* simple linear scaling -- allows more control */
129 if (simtemp->eSimTempScale == esimtempLINEAR)
131 simtemp->temperatures[i] = simtemp->simtemp_low + (simtemp->simtemp_high-simtemp->simtemp_low)*temperature_lambdas[i];
133 else if (simtemp->eSimTempScale == esimtempGEOMETRIC) /* should give roughly equal acceptance for constant heat capacity . . . */
135 simtemp->temperatures[i] = simtemp->simtemp_low * pow(simtemp->simtemp_high/simtemp->simtemp_low, (1.0*i)/(ntemps-1));
137 else if (simtemp->eSimTempScale == esimtempEXPONENTIAL)
139 simtemp->temperatures[i] = simtemp->simtemp_low + (simtemp->simtemp_high-simtemp->simtemp_low)*((exp(temperature_lambdas[i])-1)/(exp(1.0)-1));
144 sprintf(errorstr, "eSimTempScale=%d not defined", simtemp->eSimTempScale);
145 gmx_fatal(FARGS, errorstr);
152 static void _low_check(gmx_bool b, char *s, warninp_t wi)
156 warning_error(wi, s);
160 static void check_nst(const char *desc_nst, int nst,
161 const char *desc_p, int *p,
166 if (*p > 0 && *p % nst != 0)
168 /* Round up to the next multiple of nst */
169 *p = ((*p)/nst + 1)*nst;
170 sprintf(buf, "%s should be a multiple of %s, changing %s to %d\n",
171 desc_p, desc_nst, desc_p, *p);
176 static gmx_bool ir_NVE(const t_inputrec *ir)
178 return ((ir->eI == eiMD || EI_VV(ir->eI)) && ir->etc == etcNO);
181 static int lcd(int n1, int n2)
186 for (i = 2; (i <= n1 && i <= n2); i++)
188 if (n1 % i == 0 && n2 % i == 0)
197 static void process_interaction_modifier(const t_inputrec *ir, int *eintmod)
199 if (*eintmod == eintmodPOTSHIFT_VERLET)
201 if (ir->cutoff_scheme == ecutsVERLET)
203 *eintmod = eintmodPOTSHIFT;
207 *eintmod = eintmodNONE;
212 void check_ir(const char *mdparin, t_inputrec *ir, t_gromppopts *opts,
214 /* Check internal consistency */
216 /* Strange macro: first one fills the err_buf, and then one can check
217 * the condition, which will print the message and increase the error
220 #define CHECK(b) _low_check(b, err_buf, wi)
221 char err_buf[256], warn_buf[STRLEN];
227 t_lambda *fep = ir->fepvals;
228 t_expanded *expand = ir->expandedvals;
230 set_warning_line(wi, mdparin, -1);
232 /* BASIC CUT-OFF STUFF */
233 if (ir->rcoulomb < 0)
235 warning_error(wi, "rcoulomb should be >= 0");
239 warning_error(wi, "rvdw should be >= 0");
242 !(ir->cutoff_scheme == ecutsVERLET && ir->verletbuf_tol > 0))
244 warning_error(wi, "rlist should be >= 0");
247 process_interaction_modifier(ir, &ir->coulomb_modifier);
248 process_interaction_modifier(ir, &ir->vdw_modifier);
250 if (ir->cutoff_scheme == ecutsGROUP)
252 /* BASIC CUT-OFF STUFF */
253 if (ir->rlist == 0 ||
254 !((EEL_MIGHT_BE_ZERO_AT_CUTOFF(ir->coulombtype) && ir->rcoulomb > ir->rlist) ||
255 (EVDW_MIGHT_BE_ZERO_AT_CUTOFF(ir->vdwtype) && ir->rvdw > ir->rlist)))
257 /* No switched potential and/or no twin-range:
258 * we can set the long-range cut-off to the maximum of the other cut-offs.
260 ir->rlistlong = max_cutoff(ir->rlist, max_cutoff(ir->rvdw, ir->rcoulomb));
262 else if (ir->rlistlong < 0)
264 ir->rlistlong = max_cutoff(ir->rlist, max_cutoff(ir->rvdw, ir->rcoulomb));
265 sprintf(warn_buf, "rlistlong was not set, setting it to %g (no buffer)",
267 warning(wi, warn_buf);
269 if (ir->rlistlong == 0 && ir->ePBC != epbcNONE)
271 warning_error(wi, "Can not have an infinite cut-off with PBC");
273 if (ir->rlistlong > 0 && (ir->rlist == 0 || ir->rlistlong < ir->rlist))
275 warning_error(wi, "rlistlong can not be shorter than rlist");
277 if (IR_TWINRANGE(*ir) && ir->nstlist <= 0)
279 warning_error(wi, "Can not have nstlist<=0 with twin-range interactions");
283 if (ir->rlistlong == ir->rlist)
287 else if (ir->rlistlong > ir->rlist && ir->nstcalclr == 0)
289 warning_error(wi, "With different cutoffs for electrostatics and VdW, nstcalclr must be -1 or a positive number");
292 if (ir->cutoff_scheme == ecutsVERLET)
296 /* Normal Verlet type neighbor-list, currently only limited feature support */
297 if (inputrec2nboundeddim(ir) < 3)
299 warning_error(wi, "With Verlet lists only full pbc or pbc=xy with walls is supported");
301 if (ir->rcoulomb != ir->rvdw)
303 warning_error(wi, "With Verlet lists rcoulomb!=rvdw is not supported");
305 if (ir->vdwtype != evdwCUT)
307 warning_error(wi, "With Verlet lists only cut-off LJ interactions are supported");
309 if (!(ir->coulombtype == eelCUT ||
310 (EEL_RF(ir->coulombtype) && ir->coulombtype != eelRF_NEC) ||
311 EEL_PME(ir->coulombtype) || ir->coulombtype == eelEWALD))
313 warning_error(wi, "With Verlet lists only cut-off, reaction-field, PME and Ewald electrostatics are supported");
316 if (ir->nstlist <= 0)
318 warning_error(wi, "With Verlet lists nstlist should be larger than 0");
321 if (ir->nstlist < 10)
323 warning_note(wi, "With Verlet lists the optimal nstlist is >= 10, with GPUs >= 20. Note that with the Verlet scheme, nstlist has no effect on the accuracy of your simulation.");
326 rc_max = max(ir->rvdw, ir->rcoulomb);
328 if (ir->verletbuf_tol <= 0)
330 if (ir->verletbuf_tol == 0)
332 warning_error(wi, "Can not have Verlet buffer tolerance of exactly 0");
335 if (ir->rlist < rc_max)
337 warning_error(wi, "With verlet lists rlist can not be smaller than rvdw or rcoulomb");
340 if (ir->rlist == rc_max && ir->nstlist > 1)
342 warning_note(wi, "rlist is equal to rvdw and/or rcoulomb: there is no explicit Verlet buffer. The cluster pair list does have a buffering effect, but choosing a larger rlist might be necessary for good energy conservation.");
347 if (ir->rlist > rc_max)
349 warning_note(wi, "You have set rlist larger than the interaction cut-off, but you also have verlet-buffer-tolerance > 0. Will set rlist using verlet-buffer-tolerance.");
352 if (ir->nstlist == 1)
354 /* No buffer required */
359 if (EI_DYNAMICS(ir->eI))
361 if (EI_MD(ir->eI) && ir->etc == etcNO)
363 warning_error(wi, "Temperature coupling is required for calculating rlist using the energy tolerance with verlet-buffer-tolerance > 0. Either use temperature coupling or set rlist yourself together with verlet-buffer-tolerance = -1.");
366 if (inputrec2nboundeddim(ir) < 3)
368 warning_error(wi, "The box volume is required for calculating rlist from the energy drift with verlet-buffer-tolerance > 0. You are using at least one unbounded dimension, so no volume can be computed. Either use a finite box, or set rlist yourself together with verlet-buffer-tolerance = -1.");
370 /* Set rlist temporarily so we can continue processing */
375 /* Set the buffer to 5% of the cut-off */
376 ir->rlist = 1.05*rc_max;
381 /* No twin-range calculations with Verlet lists */
382 ir->rlistlong = ir->rlist;
385 if (ir->nstcalclr == -1)
387 /* if rlist=rlistlong, this will later be changed to nstcalclr=0 */
388 ir->nstcalclr = ir->nstlist;
390 else if (ir->nstcalclr > 0)
392 if (ir->nstlist > 0 && (ir->nstlist % ir->nstcalclr != 0))
394 warning_error(wi, "nstlist must be evenly divisible by nstcalclr. Use nstcalclr = -1 to automatically follow nstlist");
397 else if (ir->nstcalclr < -1)
399 warning_error(wi, "nstcalclr must be a positive number (divisor of nstcalclr), or -1 to follow nstlist.");
402 if (EEL_PME(ir->coulombtype) && ir->rcoulomb > ir->rvdw && ir->nstcalclr > 1)
404 warning_error(wi, "When used with PME, the long-range component of twin-range interactions must be updated every step (nstcalclr)");
407 /* GENERAL INTEGRATOR STUFF */
408 if (!(ir->eI == eiMD || EI_VV(ir->eI)))
412 if (ir->eI == eiVVAK)
414 sprintf(warn_buf, "Integrator method %s is implemented primarily for validation purposes; for molecular dynamics, you should probably be using %s or %s", ei_names[eiVVAK], ei_names[eiMD], ei_names[eiVV]);
415 warning_note(wi, warn_buf);
417 if (!EI_DYNAMICS(ir->eI))
421 if (EI_DYNAMICS(ir->eI))
423 if (ir->nstcalcenergy < 0)
425 ir->nstcalcenergy = ir_optimal_nstcalcenergy(ir);
426 if (ir->nstenergy != 0 && ir->nstenergy < ir->nstcalcenergy)
428 /* nstcalcenergy larger than nstener does not make sense.
429 * We ideally want nstcalcenergy=nstener.
433 ir->nstcalcenergy = lcd(ir->nstenergy, ir->nstlist);
437 ir->nstcalcenergy = ir->nstenergy;
441 else if ( (ir->nstenergy > 0 && ir->nstcalcenergy > ir->nstenergy) ||
442 (ir->efep != efepNO && ir->fepvals->nstdhdl > 0 &&
443 (ir->nstcalcenergy > ir->fepvals->nstdhdl) ) )
446 const char *nsten = "nstenergy";
447 const char *nstdh = "nstdhdl";
448 const char *min_name = nsten;
449 int min_nst = ir->nstenergy;
451 /* find the smallest of ( nstenergy, nstdhdl ) */
452 if (ir->efep != efepNO && ir->fepvals->nstdhdl > 0 &&
453 (ir->fepvals->nstdhdl < ir->nstenergy) )
455 min_nst = ir->fepvals->nstdhdl;
458 /* If the user sets nstenergy small, we should respect that */
460 "Setting nstcalcenergy (%d) equal to %s (%d)",
461 ir->nstcalcenergy, min_name, min_nst);
462 warning_note(wi, warn_buf);
463 ir->nstcalcenergy = min_nst;
466 if (ir->epc != epcNO)
468 if (ir->nstpcouple < 0)
470 ir->nstpcouple = ir_optimal_nstpcouple(ir);
473 if (IR_TWINRANGE(*ir))
475 check_nst("nstlist", ir->nstlist,
476 "nstcalcenergy", &ir->nstcalcenergy, wi);
477 if (ir->epc != epcNO)
479 check_nst("nstlist", ir->nstlist,
480 "nstpcouple", &ir->nstpcouple, wi);
484 if (ir->nstcalcenergy > 0)
486 if (ir->efep != efepNO)
488 /* nstdhdl should be a multiple of nstcalcenergy */
489 check_nst("nstcalcenergy", ir->nstcalcenergy,
490 "nstdhdl", &ir->fepvals->nstdhdl, wi);
491 /* nstexpanded should be a multiple of nstcalcenergy */
492 check_nst("nstcalcenergy", ir->nstcalcenergy,
493 "nstexpanded", &ir->expandedvals->nstexpanded, wi);
495 /* for storing exact averages nstenergy should be
496 * a multiple of nstcalcenergy
498 check_nst("nstcalcenergy", ir->nstcalcenergy,
499 "nstenergy", &ir->nstenergy, wi);
504 if ((EI_SD(ir->eI) || ir->eI == eiBD) &&
505 ir->bContinuation && ir->ld_seed != -1)
507 warning_note(wi, "You are doing a continuation with SD or BD, make sure that ld_seed is different from the previous run (using ld_seed=-1 will ensure this)");
513 sprintf(err_buf, "TPI only works with pbc = %s", epbc_names[epbcXYZ]);
514 CHECK(ir->ePBC != epbcXYZ);
515 sprintf(err_buf, "TPI only works with ns = %s", ens_names[ensGRID]);
516 CHECK(ir->ns_type != ensGRID);
517 sprintf(err_buf, "with TPI nstlist should be larger than zero");
518 CHECK(ir->nstlist <= 0);
519 sprintf(err_buf, "TPI does not work with full electrostatics other than PME");
520 CHECK(EEL_FULL(ir->coulombtype) && !EEL_PME(ir->coulombtype));
524 if ( (opts->nshake > 0) && (opts->bMorse) )
527 "Using morse bond-potentials while constraining bonds is useless");
528 warning(wi, warn_buf);
531 if ((EI_SD(ir->eI) || ir->eI == eiBD) &&
532 ir->bContinuation && ir->ld_seed != -1)
534 warning_note(wi, "You are doing a continuation with SD or BD, make sure that ld_seed is different from the previous run (using ld_seed=-1 will ensure this)");
536 /* verify simulated tempering options */
540 gmx_bool bAllTempZero = TRUE;
541 for (i = 0; i < fep->n_lambda; i++)
543 sprintf(err_buf, "Entry %d for %s must be between 0 and 1, instead is %g", i, efpt_names[efptTEMPERATURE], fep->all_lambda[efptTEMPERATURE][i]);
544 CHECK((fep->all_lambda[efptTEMPERATURE][i] < 0) || (fep->all_lambda[efptTEMPERATURE][i] > 1));
545 if (fep->all_lambda[efptTEMPERATURE][i] > 0)
547 bAllTempZero = FALSE;
550 sprintf(err_buf, "if simulated tempering is on, temperature-lambdas may not be all zero");
551 CHECK(bAllTempZero == TRUE);
553 sprintf(err_buf, "Simulated tempering is currently only compatible with md-vv");
554 CHECK(ir->eI != eiVV);
556 /* check compatability of the temperature coupling with simulated tempering */
558 if (ir->etc == etcNOSEHOOVER)
560 sprintf(warn_buf, "Nose-Hoover based temperature control such as [%s] my not be entirelyconsistent with simulated tempering", etcoupl_names[ir->etc]);
561 warning_note(wi, warn_buf);
564 /* check that the temperatures make sense */
566 sprintf(err_buf, "Higher simulated tempering temperature (%g) must be >= than the simulated tempering lower temperature (%g)", ir->simtempvals->simtemp_high, ir->simtempvals->simtemp_low);
567 CHECK(ir->simtempvals->simtemp_high <= ir->simtempvals->simtemp_low);
569 sprintf(err_buf, "Higher simulated tempering temperature (%g) must be >= zero", ir->simtempvals->simtemp_high);
570 CHECK(ir->simtempvals->simtemp_high <= 0);
572 sprintf(err_buf, "Lower simulated tempering temperature (%g) must be >= zero", ir->simtempvals->simtemp_low);
573 CHECK(ir->simtempvals->simtemp_low <= 0);
576 /* verify free energy options */
578 if (ir->efep != efepNO)
581 sprintf(err_buf, "The soft-core power is %d and can only be 1 or 2",
583 CHECK(fep->sc_alpha != 0 && fep->sc_power != 1 && fep->sc_power != 2);
585 sprintf(err_buf, "The soft-core sc-r-power is %d and can only be 6 or 48",
586 (int)fep->sc_r_power);
587 CHECK(fep->sc_alpha != 0 && fep->sc_r_power != 6.0 && fep->sc_r_power != 48.0);
589 sprintf(err_buf, "Can't use postive delta-lambda (%g) if initial state/lambda does not start at zero", fep->delta_lambda);
590 CHECK(fep->delta_lambda > 0 && ((fep->init_fep_state > 0) || (fep->init_lambda > 0)));
592 sprintf(err_buf, "Can't use postive delta-lambda (%g) with expanded ensemble simulations", fep->delta_lambda);
593 CHECK(fep->delta_lambda > 0 && (ir->efep == efepEXPANDED));
595 sprintf(err_buf, "Can only use expanded ensemble with md-vv for now; should be supported for other integrators in 5.0");
596 CHECK(!(EI_VV(ir->eI)) && (ir->efep == efepEXPANDED));
598 sprintf(err_buf, "Free-energy not implemented for Ewald");
599 CHECK(ir->coulombtype == eelEWALD);
601 /* check validty of lambda inputs */
602 if (fep->n_lambda == 0)
604 /* Clear output in case of no states:*/
605 sprintf(err_buf, "init-lambda-state set to %d: no lambda states are defined.", fep->init_fep_state);
606 CHECK((fep->init_fep_state >= 0) && (fep->n_lambda == 0));
610 sprintf(err_buf, "initial thermodynamic state %d does not exist, only goes to %d", fep->init_fep_state, fep->n_lambda-1);
611 CHECK((fep->init_fep_state >= fep->n_lambda));
614 sprintf(err_buf, "Lambda state must be set, either with init-lambda-state or with init-lambda");
615 CHECK((fep->init_fep_state < 0) && (fep->init_lambda < 0));
617 sprintf(err_buf, "init-lambda=%g while init-lambda-state=%d. Lambda state must be set either with init-lambda-state or with init-lambda, but not both",
618 fep->init_lambda, fep->init_fep_state);
619 CHECK((fep->init_fep_state >= 0) && (fep->init_lambda >= 0));
623 if ((fep->init_lambda >= 0) && (fep->delta_lambda == 0))
627 for (i = 0; i < efptNR; i++)
629 if (fep->separate_dvdl[i])
634 if (n_lambda_terms > 1)
636 sprintf(warn_buf, "If lambda vector states (fep-lambdas, coul-lambdas etc.) are set, don't use init-lambda to set lambda state (except for slow growth). Use init-lambda-state instead.");
637 warning(wi, warn_buf);
640 if (n_lambda_terms < 2 && fep->n_lambda > 0)
643 "init-lambda is deprecated for setting lambda state (except for slow growth). Use init-lambda-state instead.");
647 for (j = 0; j < efptNR; j++)
649 for (i = 0; i < fep->n_lambda; i++)
651 sprintf(err_buf, "Entry %d for %s must be between 0 and 1, instead is %g", i, efpt_names[j], fep->all_lambda[j][i]);
652 CHECK((fep->all_lambda[j][i] < 0) || (fep->all_lambda[j][i] > 1));
656 if ((fep->sc_alpha > 0) && (!fep->bScCoul))
658 for (i = 0; i < fep->n_lambda; i++)
660 sprintf(err_buf, "For state %d, vdw-lambdas (%f) is changing with vdw softcore, while coul-lambdas (%f) is nonzero without coulomb softcore: this will lead to crashes, and is not supported.", i, fep->all_lambda[efptVDW][i],
661 fep->all_lambda[efptCOUL][i]);
662 CHECK((fep->sc_alpha > 0) &&
663 (((fep->all_lambda[efptCOUL][i] > 0.0) &&
664 (fep->all_lambda[efptCOUL][i] < 1.0)) &&
665 ((fep->all_lambda[efptVDW][i] > 0.0) &&
666 (fep->all_lambda[efptVDW][i] < 1.0))));
670 if ((fep->bScCoul) && (EEL_PME(ir->coulombtype)))
672 real sigma, lambda, r_sc;
675 /* Maximum estimate for A and B charges equal with lambda power 1 */
677 r_sc = pow(lambda*fep->sc_alpha*pow(sigma/ir->rcoulomb, fep->sc_r_power) + 1.0, 1.0/fep->sc_r_power);
678 sprintf(warn_buf, "With PME there is a minor soft core effect present at the cut-off, proportional to (LJsigma/rcoulomb)^%g. This could have a minor effect on energy conservation, but usually other effects dominate. With a common sigma value of %g nm the fraction of the particle-particle potential at the cut-off at lambda=%g is around %.1e, while ewald-rtol is %.1e.",
680 sigma, lambda, r_sc - 1.0, ir->ewald_rtol);
681 warning_note(wi, warn_buf);
684 /* Free Energy Checks -- In an ideal world, slow growth and FEP would
685 be treated differently, but that's the next step */
687 for (i = 0; i < efptNR; i++)
689 for (j = 0; j < fep->n_lambda; j++)
691 sprintf(err_buf, "%s[%d] must be between 0 and 1", efpt_names[i], j);
692 CHECK((fep->all_lambda[i][j] < 0) || (fep->all_lambda[i][j] > 1));
697 if ((ir->bSimTemp) || (ir->efep == efepEXPANDED))
700 expand = ir->expandedvals;
702 /* checking equilibration of weights inputs for validity */
704 sprintf(err_buf, "weight-equil-number-all-lambda (%d) is ignored if lmc-weights-equil is not equal to %s",
705 expand->equil_n_at_lam, elmceq_names[elmceqNUMATLAM]);
706 CHECK((expand->equil_n_at_lam > 0) && (expand->elmceq != elmceqNUMATLAM));
708 sprintf(err_buf, "weight-equil-number-samples (%d) is ignored if lmc-weights-equil is not equal to %s",
709 expand->equil_samples, elmceq_names[elmceqSAMPLES]);
710 CHECK((expand->equil_samples > 0) && (expand->elmceq != elmceqSAMPLES));
712 sprintf(err_buf, "weight-equil-number-steps (%d) is ignored if lmc-weights-equil is not equal to %s",
713 expand->equil_steps, elmceq_names[elmceqSTEPS]);
714 CHECK((expand->equil_steps > 0) && (expand->elmceq != elmceqSTEPS));
716 sprintf(err_buf, "weight-equil-wl-delta (%d) is ignored if lmc-weights-equil is not equal to %s",
717 expand->equil_samples, elmceq_names[elmceqWLDELTA]);
718 CHECK((expand->equil_wl_delta > 0) && (expand->elmceq != elmceqWLDELTA));
720 sprintf(err_buf, "weight-equil-count-ratio (%f) is ignored if lmc-weights-equil is not equal to %s",
721 expand->equil_ratio, elmceq_names[elmceqRATIO]);
722 CHECK((expand->equil_ratio > 0) && (expand->elmceq != elmceqRATIO));
724 sprintf(err_buf, "weight-equil-number-all-lambda (%d) must be a positive integer if lmc-weights-equil=%s",
725 expand->equil_n_at_lam, elmceq_names[elmceqNUMATLAM]);
726 CHECK((expand->equil_n_at_lam <= 0) && (expand->elmceq == elmceqNUMATLAM));
728 sprintf(err_buf, "weight-equil-number-samples (%d) must be a positive integer if lmc-weights-equil=%s",
729 expand->equil_samples, elmceq_names[elmceqSAMPLES]);
730 CHECK((expand->equil_samples <= 0) && (expand->elmceq == elmceqSAMPLES));
732 sprintf(err_buf, "weight-equil-number-steps (%d) must be a positive integer if lmc-weights-equil=%s",
733 expand->equil_steps, elmceq_names[elmceqSTEPS]);
734 CHECK((expand->equil_steps <= 0) && (expand->elmceq == elmceqSTEPS));
736 sprintf(err_buf, "weight-equil-wl-delta (%f) must be > 0 if lmc-weights-equil=%s",
737 expand->equil_wl_delta, elmceq_names[elmceqWLDELTA]);
738 CHECK((expand->equil_wl_delta <= 0) && (expand->elmceq == elmceqWLDELTA));
740 sprintf(err_buf, "weight-equil-count-ratio (%f) must be > 0 if lmc-weights-equil=%s",
741 expand->equil_ratio, elmceq_names[elmceqRATIO]);
742 CHECK((expand->equil_ratio <= 0) && (expand->elmceq == elmceqRATIO));
744 sprintf(err_buf, "lmc-weights-equil=%s only possible when lmc-stats = %s or lmc-stats %s",
745 elmceq_names[elmceqWLDELTA], elamstats_names[elamstatsWL], elamstats_names[elamstatsWWL]);
746 CHECK((expand->elmceq == elmceqWLDELTA) && (!EWL(expand->elamstats)));
748 sprintf(err_buf, "lmc-repeats (%d) must be greater than 0", expand->lmc_repeats);
749 CHECK((expand->lmc_repeats <= 0));
750 sprintf(err_buf, "minimum-var-min (%d) must be greater than 0", expand->minvarmin);
751 CHECK((expand->minvarmin <= 0));
752 sprintf(err_buf, "weight-c-range (%d) must be greater or equal to 0", expand->c_range);
753 CHECK((expand->c_range < 0));
754 sprintf(err_buf, "init-lambda-state (%d) must be zero if lmc-forced-nstart (%d)> 0 and lmc-move != 'no'",
755 fep->init_fep_state, expand->lmc_forced_nstart);
756 CHECK((fep->init_fep_state != 0) && (expand->lmc_forced_nstart > 0) && (expand->elmcmove != elmcmoveNO));
757 sprintf(err_buf, "lmc-forced-nstart (%d) must not be negative", expand->lmc_forced_nstart);
758 CHECK((expand->lmc_forced_nstart < 0));
759 sprintf(err_buf, "init-lambda-state (%d) must be in the interval [0,number of lambdas)", fep->init_fep_state);
760 CHECK((fep->init_fep_state < 0) || (fep->init_fep_state >= fep->n_lambda));
762 sprintf(err_buf, "init-wl-delta (%f) must be greater than or equal to 0", expand->init_wl_delta);
763 CHECK((expand->init_wl_delta < 0));
764 sprintf(err_buf, "wl-ratio (%f) must be between 0 and 1", expand->wl_ratio);
765 CHECK((expand->wl_ratio <= 0) || (expand->wl_ratio >= 1));
766 sprintf(err_buf, "wl-scale (%f) must be between 0 and 1", expand->wl_scale);
767 CHECK((expand->wl_scale <= 0) || (expand->wl_scale >= 1));
769 /* if there is no temperature control, we need to specify an MC temperature */
770 sprintf(err_buf, "If there is no temperature control, and lmc-mcmove!= 'no',mc_temperature must be set to a positive number");
771 if (expand->nstTij > 0)
773 sprintf(err_buf, "nst-transition-matrix (%d) must be an integer multiple of nstlog (%d)",
774 expand->nstTij, ir->nstlog);
775 CHECK((mod(expand->nstTij, ir->nstlog) != 0));
780 sprintf(err_buf, "walls only work with pbc=%s", epbc_names[epbcXY]);
781 CHECK(ir->nwall && ir->ePBC != epbcXY);
784 if (ir->ePBC != epbcXYZ && ir->nwall != 2)
786 if (ir->ePBC == epbcNONE)
788 if (ir->epc != epcNO)
790 warning(wi, "Turning off pressure coupling for vacuum system");
796 sprintf(err_buf, "Can not have pressure coupling with pbc=%s",
797 epbc_names[ir->ePBC]);
798 CHECK(ir->epc != epcNO);
800 sprintf(err_buf, "Can not have Ewald with pbc=%s", epbc_names[ir->ePBC]);
801 CHECK(EEL_FULL(ir->coulombtype));
803 sprintf(err_buf, "Can not have dispersion correction with pbc=%s",
804 epbc_names[ir->ePBC]);
805 CHECK(ir->eDispCorr != edispcNO);
808 if (ir->rlist == 0.0)
810 sprintf(err_buf, "can only have neighborlist cut-off zero (=infinite)\n"
811 "with coulombtype = %s or coulombtype = %s\n"
812 "without periodic boundary conditions (pbc = %s) and\n"
813 "rcoulomb and rvdw set to zero",
814 eel_names[eelCUT], eel_names[eelUSER], epbc_names[epbcNONE]);
815 CHECK(((ir->coulombtype != eelCUT) && (ir->coulombtype != eelUSER)) ||
816 (ir->ePBC != epbcNONE) ||
817 (ir->rcoulomb != 0.0) || (ir->rvdw != 0.0));
821 warning_error(wi, "Can not have heuristic neighborlist updates without cut-off");
825 warning_note(wi, "Simulating without cut-offs is usually (slightly) faster with nstlist=0, nstype=simple and particle decomposition");
830 if (ir->nstcomm == 0)
832 ir->comm_mode = ecmNO;
834 if (ir->comm_mode != ecmNO)
838 warning(wi, "If you want to remove the rotation around the center of mass, you should set comm_mode = Angular instead of setting nstcomm < 0. nstcomm is modified to its absolute value");
839 ir->nstcomm = abs(ir->nstcomm);
842 if (ir->nstcalcenergy > 0 && ir->nstcomm < ir->nstcalcenergy)
844 warning_note(wi, "nstcomm < nstcalcenergy defeats the purpose of nstcalcenergy, setting nstcomm to nstcalcenergy");
845 ir->nstcomm = ir->nstcalcenergy;
848 if (ir->comm_mode == ecmANGULAR)
850 sprintf(err_buf, "Can not remove the rotation around the center of mass with periodic molecules");
851 CHECK(ir->bPeriodicMols);
852 if (ir->ePBC != epbcNONE)
854 warning(wi, "Removing the rotation around the center of mass in a periodic system (this is not a problem when you have only one molecule).");
859 if (EI_STATE_VELOCITY(ir->eI) && ir->ePBC == epbcNONE && ir->comm_mode != ecmANGULAR)
861 warning_note(wi, "Tumbling and or flying ice-cubes: We are not removing rotation around center of mass in a non-periodic system. You should probably set comm_mode = ANGULAR.");
864 sprintf(err_buf, "Twin-range neighbour searching (NS) with simple NS"
865 " algorithm not implemented");
866 CHECK(((ir->rcoulomb > ir->rlist) || (ir->rvdw > ir->rlist))
867 && (ir->ns_type == ensSIMPLE));
869 /* TEMPERATURE COUPLING */
870 if (ir->etc == etcYES)
872 ir->etc = etcBERENDSEN;
873 warning_note(wi, "Old option for temperature coupling given: "
874 "changing \"yes\" to \"Berendsen\"\n");
877 if ((ir->etc == etcNOSEHOOVER) || (ir->epc == epcMTTK))
879 if (ir->opts.nhchainlength < 1)
881 sprintf(warn_buf, "number of Nose-Hoover chains (currently %d) cannot be less than 1,reset to 1\n", ir->opts.nhchainlength);
882 ir->opts.nhchainlength = 1;
883 warning(wi, warn_buf);
886 if (ir->etc == etcNOSEHOOVER && !EI_VV(ir->eI) && ir->opts.nhchainlength > 1)
888 warning_note(wi, "leapfrog does not yet support Nose-Hoover chains, nhchainlength reset to 1");
889 ir->opts.nhchainlength = 1;
894 ir->opts.nhchainlength = 0;
897 if (ir->eI == eiVVAK)
899 sprintf(err_buf, "%s implemented primarily for validation, and requires nsttcouple = 1 and nstpcouple = 1.",
901 CHECK((ir->nsttcouple != 1) || (ir->nstpcouple != 1));
904 if (ETC_ANDERSEN(ir->etc))
906 sprintf(err_buf, "%s temperature control not supported for integrator %s.", etcoupl_names[ir->etc], ei_names[ir->eI]);
907 CHECK(!(EI_VV(ir->eI)));
909 for (i = 0; i < ir->opts.ngtc; i++)
911 sprintf(err_buf, "all tau_t must currently be equal using Andersen temperature control, violated for group %d", i);
912 CHECK(ir->opts.tau_t[0] != ir->opts.tau_t[i]);
913 sprintf(err_buf, "all tau_t must be postive using Andersen temperature control, tau_t[%d]=%10.6f",
914 i, ir->opts.tau_t[i]);
915 CHECK(ir->opts.tau_t[i] < 0);
917 if (ir->nstcomm > 0 && (ir->etc == etcANDERSEN))
919 sprintf(warn_buf, "Center of mass removal not necessary for %s. All velocities of coupled groups are rerandomized periodically, so flying ice cube errors will not occur.", etcoupl_names[ir->etc]);
920 warning_note(wi, warn_buf);
923 sprintf(err_buf, "nstcomm must be 1, not %d for %s, as velocities of atoms in coupled groups are randomized every time step", ir->nstcomm, etcoupl_names[ir->etc]);
924 CHECK(ir->nstcomm > 1 && (ir->etc == etcANDERSEN));
926 for (i = 0; i < ir->opts.ngtc; i++)
928 int nsteps = (int)(ir->opts.tau_t[i]/ir->delta_t);
929 sprintf(err_buf, "tau_t/delta_t for group %d for temperature control method %s must be a multiple of nstcomm (%d), as velocities of atoms in coupled groups are randomized every time step. The input tau_t (%8.3f) leads to %d steps per randomization", i, etcoupl_names[ir->etc], ir->nstcomm, ir->opts.tau_t[i], nsteps);
930 CHECK((nsteps % ir->nstcomm) && (ir->etc == etcANDERSENMASSIVE));
933 if (ir->etc == etcBERENDSEN)
935 sprintf(warn_buf, "The %s thermostat does not generate the correct kinetic energy distribution. You might want to consider using the %s thermostat.",
936 ETCOUPLTYPE(ir->etc), ETCOUPLTYPE(etcVRESCALE));
937 warning_note(wi, warn_buf);
940 if ((ir->etc == etcNOSEHOOVER || ETC_ANDERSEN(ir->etc))
941 && ir->epc == epcBERENDSEN)
943 sprintf(warn_buf, "Using Berendsen pressure coupling invalidates the "
944 "true ensemble for the thermostat");
945 warning(wi, warn_buf);
948 /* PRESSURE COUPLING */
949 if (ir->epc == epcISOTROPIC)
951 ir->epc = epcBERENDSEN;
952 warning_note(wi, "Old option for pressure coupling given: "
953 "changing \"Isotropic\" to \"Berendsen\"\n");
956 if (ir->epc != epcNO)
958 dt_pcoupl = ir->nstpcouple*ir->delta_t;
960 sprintf(err_buf, "tau-p must be > 0 instead of %g\n", ir->tau_p);
961 CHECK(ir->tau_p <= 0);
963 if (ir->tau_p/dt_pcoupl < pcouple_min_integration_steps(ir->epc))
965 sprintf(warn_buf, "For proper integration of the %s barostat, tau-p (%g) should be at least %d times larger than nstpcouple*dt (%g)",
966 EPCOUPLTYPE(ir->epc), ir->tau_p, pcouple_min_integration_steps(ir->epc), dt_pcoupl);
967 warning(wi, warn_buf);
970 sprintf(err_buf, "compressibility must be > 0 when using pressure"
971 " coupling %s\n", EPCOUPLTYPE(ir->epc));
972 CHECK(ir->compress[XX][XX] < 0 || ir->compress[YY][YY] < 0 ||
973 ir->compress[ZZ][ZZ] < 0 ||
974 (trace(ir->compress) == 0 && ir->compress[YY][XX] <= 0 &&
975 ir->compress[ZZ][XX] <= 0 && ir->compress[ZZ][YY] <= 0));
977 if (epcPARRINELLORAHMAN == ir->epc && opts->bGenVel)
980 "You are generating velocities so I am assuming you "
981 "are equilibrating a system. You are using "
982 "%s pressure coupling, but this can be "
983 "unstable for equilibration. If your system crashes, try "
984 "equilibrating first with Berendsen pressure coupling. If "
985 "you are not equilibrating the system, you can probably "
986 "ignore this warning.",
987 epcoupl_names[ir->epc]);
988 warning(wi, warn_buf);
996 if ((ir->epc != epcBERENDSEN) && (ir->epc != epcMTTK))
998 warning_error(wi, "for md-vv and md-vv-avek, can only use Berendsen and Martyna-Tuckerman-Tobias-Klein (MTTK) equations for pressure control; MTTK is equivalent to Parrinello-Rahman.");
1003 /* ELECTROSTATICS */
1004 /* More checks are in triple check (grompp.c) */
1006 if (ir->coulombtype == eelSWITCH)
1008 sprintf(warn_buf, "coulombtype = %s is only for testing purposes and can lead to serious "
1009 "artifacts, advice: use coulombtype = %s",
1010 eel_names[ir->coulombtype],
1011 eel_names[eelRF_ZERO]);
1012 warning(wi, warn_buf);
1015 if (ir->epsilon_r != 1 && ir->implicit_solvent == eisGBSA)
1017 sprintf(warn_buf, "epsilon-r = %g with GB implicit solvent, will use this value for inner dielectric", ir->epsilon_r);
1018 warning_note(wi, warn_buf);
1021 if (EEL_RF(ir->coulombtype) && ir->epsilon_rf == 1 && ir->epsilon_r != 1)
1023 sprintf(warn_buf, "epsilon-r = %g and epsilon-rf = 1 with reaction field, proceeding assuming old format and exchanging epsilon-r and epsilon-rf", ir->epsilon_r);
1024 warning(wi, warn_buf);
1025 ir->epsilon_rf = ir->epsilon_r;
1026 ir->epsilon_r = 1.0;
1029 if (getenv("GALACTIC_DYNAMICS") == NULL)
1031 sprintf(err_buf, "epsilon-r must be >= 0 instead of %g\n", ir->epsilon_r);
1032 CHECK(ir->epsilon_r < 0);
1035 if (EEL_RF(ir->coulombtype))
1037 /* reaction field (at the cut-off) */
1039 if (ir->coulombtype == eelRF_ZERO)
1041 sprintf(warn_buf, "With coulombtype = %s, epsilon-rf must be 0, assuming you meant epsilon_rf=0",
1042 eel_names[ir->coulombtype]);
1043 CHECK(ir->epsilon_rf != 0);
1044 ir->epsilon_rf = 0.0;
1047 sprintf(err_buf, "epsilon-rf must be >= epsilon-r");
1048 CHECK((ir->epsilon_rf < ir->epsilon_r && ir->epsilon_rf != 0) ||
1049 (ir->epsilon_r == 0));
1050 if (ir->epsilon_rf == ir->epsilon_r)
1052 sprintf(warn_buf, "Using epsilon-rf = epsilon-r with %s does not make sense",
1053 eel_names[ir->coulombtype]);
1054 warning(wi, warn_buf);
1057 /* Allow rlist>rcoulomb for tabulated long range stuff. This just
1058 * means the interaction is zero outside rcoulomb, but it helps to
1059 * provide accurate energy conservation.
1061 if (EEL_MIGHT_BE_ZERO_AT_CUTOFF(ir->coulombtype))
1063 if (EEL_SWITCHED(ir->coulombtype))
1066 "With coulombtype = %s rcoulomb_switch must be < rcoulomb. Or, better: Use the potential modifier options!",
1067 eel_names[ir->coulombtype]);
1068 CHECK(ir->rcoulomb_switch >= ir->rcoulomb);
1071 else if (ir->coulombtype == eelCUT || EEL_RF(ir->coulombtype))
1073 if (ir->cutoff_scheme == ecutsGROUP && ir->coulomb_modifier == eintmodNONE)
1075 sprintf(err_buf, "With coulombtype = %s, rcoulomb should be >= rlist unless you use a potential modifier",
1076 eel_names[ir->coulombtype]);
1077 CHECK(ir->rlist > ir->rcoulomb);
1081 if (ir->coulombtype == eelSWITCH || ir->coulombtype == eelSHIFT ||
1082 ir->vdwtype == evdwSWITCH || ir->vdwtype == evdwSHIFT)
1085 "The switch/shift interaction settings are just for compatibility; you will get better "
1086 "performance from applying potential modifiers to your interactions!\n");
1087 warning_note(wi, warn_buf);
1090 if (ir->coulombtype == eelPMESWITCH)
1092 if (ir->rcoulomb_switch/ir->rcoulomb < 0.9499)
1094 sprintf(warn_buf, "The switching range for %s should be 5%% or less, energy conservation will be good anyhow, since ewald_rtol = %g",
1095 eel_names[ir->coulombtype],
1097 warning(wi, warn_buf);
1101 if (EEL_FULL(ir->coulombtype))
1103 if (ir->coulombtype == eelPMESWITCH || ir->coulombtype == eelPMEUSER ||
1104 ir->coulombtype == eelPMEUSERSWITCH)
1106 sprintf(err_buf, "With coulombtype = %s, rcoulomb must be <= rlist",
1107 eel_names[ir->coulombtype]);
1108 CHECK(ir->rcoulomb > ir->rlist);
1110 else if (ir->cutoff_scheme == ecutsGROUP && ir->coulomb_modifier == eintmodNONE)
1112 if (ir->coulombtype == eelPME || ir->coulombtype == eelP3M_AD)
1115 "With coulombtype = %s (without modifier), rcoulomb must be equal to rlist,\n"
1116 "or rlistlong if nstcalclr=1. For optimal energy conservation,consider using\n"
1117 "a potential modifier.", eel_names[ir->coulombtype]);
1118 if (ir->nstcalclr == 1)
1120 CHECK(ir->rcoulomb != ir->rlist && ir->rcoulomb != ir->rlistlong);
1124 CHECK(ir->rcoulomb != ir->rlist);
1130 if (EVDW_PME(ir->vdwtype))
1132 if (EVDW_MIGHT_BE_ZERO_AT_CUTOFF(ir->vdwtype))
1134 sprintf(err_buf, "With vdwtype = %s, rvdw must be <= rlist",
1135 evdw_names[ir->vdwtype]);
1136 CHECK(ir->rvdw > ir->rlist);
1141 "With vdwtype = %s, rvdw must be equal to rlist\n",
1142 evdw_names[ir->vdwtype]);
1143 CHECK(ir->rvdw != ir->rlist);
1147 if (EEL_PME(ir->coulombtype) || EVDW_PME(ir->vdwtype))
1149 if (ir->pme_order < 3)
1151 warning_error(wi, "pme-order can not be smaller than 3");
1155 if (ir->nwall == 2 && EEL_FULL(ir->coulombtype))
1157 if (ir->ewald_geometry == eewg3D)
1159 sprintf(warn_buf, "With pbc=%s you should use ewald-geometry=%s",
1160 epbc_names[ir->ePBC], eewg_names[eewg3DC]);
1161 warning(wi, warn_buf);
1163 /* This check avoids extra pbc coding for exclusion corrections */
1164 sprintf(err_buf, "wall-ewald-zfac should be >= 2");
1165 CHECK(ir->wall_ewald_zfac < 2);
1168 if (EVDW_SWITCHED(ir->vdwtype))
1170 sprintf(err_buf, "With vdwtype = %s rvdw-switch must be < rvdw. Or, better - use a potential modifier.",
1171 evdw_names[ir->vdwtype]);
1172 CHECK(ir->rvdw_switch >= ir->rvdw);
1174 else if (ir->vdwtype == evdwCUT)
1176 if (ir->cutoff_scheme == ecutsGROUP && ir->vdw_modifier == eintmodNONE)
1178 sprintf(err_buf, "With vdwtype = %s, rvdw must be >= rlist unless you use a potential modifier", evdw_names[ir->vdwtype]);
1179 CHECK(ir->rlist > ir->rvdw);
1182 if (ir->cutoff_scheme == ecutsGROUP)
1184 if (((ir->coulomb_modifier != eintmodNONE && ir->rcoulomb == ir->rlist) ||
1185 (ir->vdw_modifier != eintmodNONE && ir->rvdw == ir->rlist)) &&
1188 warning_note(wi, "With exact cut-offs, rlist should be "
1189 "larger than rcoulomb and rvdw, so that there "
1190 "is a buffer region for particle motion "
1191 "between neighborsearch steps");
1194 if (EEL_IS_ZERO_AT_CUTOFF(ir->coulombtype)
1195 && (ir->rlistlong <= ir->rcoulomb))
1197 sprintf(warn_buf, "For energy conservation with switch/shift potentials, %s should be 0.1 to 0.3 nm larger than rcoulomb.",
1198 IR_TWINRANGE(*ir) ? "rlistlong" : "rlist");
1199 warning_note(wi, warn_buf);
1201 if (EVDW_SWITCHED(ir->vdwtype) && (ir->rlistlong <= ir->rvdw))
1203 sprintf(warn_buf, "For energy conservation with switch/shift potentials, %s should be 0.1 to 0.3 nm larger than rvdw.",
1204 IR_TWINRANGE(*ir) ? "rlistlong" : "rlist");
1205 warning_note(wi, warn_buf);
1209 if (ir->vdwtype == evdwUSER && ir->eDispCorr != edispcNO)
1211 warning_note(wi, "You have selected user tables with dispersion correction, the dispersion will be corrected to -C6/r^6 beyond rvdw_switch (the tabulated interaction between rvdw_switch and rvdw will not be double counted). Make sure that you really want dispersion correction to -C6/r^6.");
1214 if (ir->nstlist == -1)
1216 sprintf(err_buf, "With nstlist=-1 rvdw and rcoulomb should be smaller than rlist to account for diffusion and possibly charge-group radii");
1217 CHECK(ir->rvdw >= ir->rlist || ir->rcoulomb >= ir->rlist);
1219 sprintf(err_buf, "nstlist can not be smaller than -1");
1220 CHECK(ir->nstlist < -1);
1222 if (ir->eI == eiLBFGS && (ir->coulombtype == eelCUT || ir->vdwtype == evdwCUT)
1225 warning(wi, "For efficient BFGS minimization, use switch/shift/pme instead of cut-off.");
1228 if (ir->eI == eiLBFGS && ir->nbfgscorr <= 0)
1230 warning(wi, "Using L-BFGS with nbfgscorr<=0 just gets you steepest descent.");
1233 /* ENERGY CONSERVATION */
1234 if (ir_NVE(ir) && ir->cutoff_scheme == ecutsGROUP)
1236 if (!EVDW_MIGHT_BE_ZERO_AT_CUTOFF(ir->vdwtype) && ir->rvdw > 0 && ir->vdw_modifier == eintmodNONE)
1238 sprintf(warn_buf, "You are using a cut-off for VdW interactions with NVE, for good energy conservation use vdwtype = %s (possibly with DispCorr)",
1239 evdw_names[evdwSHIFT]);
1240 warning_note(wi, warn_buf);
1242 if (!EEL_MIGHT_BE_ZERO_AT_CUTOFF(ir->coulombtype) && ir->rcoulomb > 0 && ir->coulomb_modifier == eintmodNONE)
1244 sprintf(warn_buf, "You are using a cut-off for electrostatics with NVE, for good energy conservation use coulombtype = %s or %s",
1245 eel_names[eelPMESWITCH], eel_names[eelRF_ZERO]);
1246 warning_note(wi, warn_buf);
1250 /* IMPLICIT SOLVENT */
1251 if (ir->coulombtype == eelGB_NOTUSED)
1253 ir->coulombtype = eelCUT;
1254 ir->implicit_solvent = eisGBSA;
1255 fprintf(stderr, "Note: Old option for generalized born electrostatics given:\n"
1256 "Changing coulombtype from \"generalized-born\" to \"cut-off\" and instead\n"
1257 "setting implicit-solvent value to \"GBSA\" in input section.\n");
1260 if (ir->sa_algorithm == esaSTILL)
1262 sprintf(err_buf, "Still SA algorithm not available yet, use %s or %s instead\n", esa_names[esaAPPROX], esa_names[esaNO]);
1263 CHECK(ir->sa_algorithm == esaSTILL);
1266 if (ir->implicit_solvent == eisGBSA)
1268 sprintf(err_buf, "With GBSA implicit solvent, rgbradii must be equal to rlist.");
1269 CHECK(ir->rgbradii != ir->rlist);
1271 if (ir->coulombtype != eelCUT)
1273 sprintf(err_buf, "With GBSA, coulombtype must be equal to %s\n", eel_names[eelCUT]);
1274 CHECK(ir->coulombtype != eelCUT);
1276 if (ir->vdwtype != evdwCUT)
1278 sprintf(err_buf, "With GBSA, vdw-type must be equal to %s\n", evdw_names[evdwCUT]);
1279 CHECK(ir->vdwtype != evdwCUT);
1281 if (ir->nstgbradii < 1)
1283 sprintf(warn_buf, "Using GBSA with nstgbradii<1, setting nstgbradii=1");
1284 warning_note(wi, warn_buf);
1287 if (ir->sa_algorithm == esaNO)
1289 sprintf(warn_buf, "No SA (non-polar) calculation requested together with GB. Are you sure this is what you want?\n");
1290 warning_note(wi, warn_buf);
1292 if (ir->sa_surface_tension < 0 && ir->sa_algorithm != esaNO)
1294 sprintf(warn_buf, "Value of sa_surface_tension is < 0. Changing it to 2.05016 or 2.25936 kJ/nm^2/mol for Still and HCT/OBC respectively\n");
1295 warning_note(wi, warn_buf);
1297 if (ir->gb_algorithm == egbSTILL)
1299 ir->sa_surface_tension = 0.0049 * CAL2JOULE * 100;
1303 ir->sa_surface_tension = 0.0054 * CAL2JOULE * 100;
1306 if (ir->sa_surface_tension == 0 && ir->sa_algorithm != esaNO)
1308 sprintf(err_buf, "Surface tension set to 0 while SA-calculation requested\n");
1309 CHECK(ir->sa_surface_tension == 0 && ir->sa_algorithm != esaNO);
1316 if (ir->cutoff_scheme != ecutsGROUP)
1318 warning_error(wi, "AdresS simulation supports only cutoff-scheme=group");
1322 warning_error(wi, "AdresS simulation supports only stochastic dynamics");
1324 if (ir->epc != epcNO)
1326 warning_error(wi, "AdresS simulation does not support pressure coupling");
1328 if (EEL_FULL(ir->coulombtype))
1330 warning_error(wi, "AdresS simulation does not support long-range electrostatics");
1335 /* count the number of text elemets separated by whitespace in a string.
1336 str = the input string
1337 maxptr = the maximum number of allowed elements
1338 ptr = the output array of pointers to the first character of each element
1339 returns: the number of elements. */
1340 int str_nelem(const char *str, int maxptr, char *ptr[])
1345 copy0 = strdup(str);
1348 while (*copy != '\0')
1352 gmx_fatal(FARGS, "Too many groups on line: '%s' (max is %d)",
1360 while ((*copy != '\0') && !isspace(*copy))
1379 /* interpret a number of doubles from a string and put them in an array,
1380 after allocating space for them.
1381 str = the input string
1382 n = the (pre-allocated) number of doubles read
1383 r = the output array of doubles. */
1384 static void parse_n_real(char *str, int *n, real **r)
1389 *n = str_nelem(str, MAXPTR, ptr);
1392 for (i = 0; i < *n; i++)
1394 (*r)[i] = strtod(ptr[i], NULL);
1398 static void do_fep_params(t_inputrec *ir, char fep_lambda[][STRLEN], char weights[STRLEN])
1401 int i, j, max_n_lambda, nweights, nfep[efptNR];
1402 t_lambda *fep = ir->fepvals;
1403 t_expanded *expand = ir->expandedvals;
1404 real **count_fep_lambdas;
1405 gmx_bool bOneLambda = TRUE;
1407 snew(count_fep_lambdas, efptNR);
1409 /* FEP input processing */
1410 /* first, identify the number of lambda values for each type.
1411 All that are nonzero must have the same number */
1413 for (i = 0; i < efptNR; i++)
1415 parse_n_real(fep_lambda[i], &(nfep[i]), &(count_fep_lambdas[i]));
1418 /* now, determine the number of components. All must be either zero, or equal. */
1421 for (i = 0; i < efptNR; i++)
1423 if (nfep[i] > max_n_lambda)
1425 max_n_lambda = nfep[i]; /* here's a nonzero one. All of them
1426 must have the same number if its not zero.*/
1431 for (i = 0; i < efptNR; i++)
1435 ir->fepvals->separate_dvdl[i] = FALSE;
1437 else if (nfep[i] == max_n_lambda)
1439 if (i != efptTEMPERATURE) /* we treat this differently -- not really a reason to compute the derivative with
1440 respect to the temperature currently */
1442 ir->fepvals->separate_dvdl[i] = TRUE;
1447 gmx_fatal(FARGS, "Number of lambdas (%d) for FEP type %s not equal to number of other types (%d)",
1448 nfep[i], efpt_names[i], max_n_lambda);
1451 /* we don't print out dhdl if the temperature is changing, since we can't correctly define dhdl in this case */
1452 ir->fepvals->separate_dvdl[efptTEMPERATURE] = FALSE;
1454 /* the number of lambdas is the number we've read in, which is either zero
1455 or the same for all */
1456 fep->n_lambda = max_n_lambda;
1458 /* allocate space for the array of lambda values */
1459 snew(fep->all_lambda, efptNR);
1460 /* if init_lambda is defined, we need to set lambda */
1461 if ((fep->init_lambda > 0) && (fep->n_lambda == 0))
1463 ir->fepvals->separate_dvdl[efptFEP] = TRUE;
1465 /* otherwise allocate the space for all of the lambdas, and transfer the data */
1466 for (i = 0; i < efptNR; i++)
1468 snew(fep->all_lambda[i], fep->n_lambda);
1469 if (nfep[i] > 0) /* if it's zero, then the count_fep_lambda arrays
1472 for (j = 0; j < fep->n_lambda; j++)
1474 fep->all_lambda[i][j] = (double)count_fep_lambdas[i][j];
1476 sfree(count_fep_lambdas[i]);
1479 sfree(count_fep_lambdas);
1481 /* "fep-vals" is either zero or the full number. If zero, we'll need to define fep-lambdas for internal
1482 bookkeeping -- for now, init_lambda */
1484 if ((nfep[efptFEP] == 0) && (fep->init_lambda >= 0))
1486 for (i = 0; i < fep->n_lambda; i++)
1488 fep->all_lambda[efptFEP][i] = fep->init_lambda;
1492 /* check to see if only a single component lambda is defined, and soft core is defined.
1493 In this case, turn on coulomb soft core */
1495 if (max_n_lambda == 0)
1501 for (i = 0; i < efptNR; i++)
1503 if ((nfep[i] != 0) && (i != efptFEP))
1509 if ((bOneLambda) && (fep->sc_alpha > 0))
1511 fep->bScCoul = TRUE;
1514 /* Fill in the others with the efptFEP if they are not explicitly
1515 specified (i.e. nfep[i] == 0). This means if fep is not defined,
1516 they are all zero. */
1518 for (i = 0; i < efptNR; i++)
1520 if ((nfep[i] == 0) && (i != efptFEP))
1522 for (j = 0; j < fep->n_lambda; j++)
1524 fep->all_lambda[i][j] = fep->all_lambda[efptFEP][j];
1530 /* make it easier if sc_r_power = 48 by increasing it to the 4th power, to be in the right scale. */
1531 if (fep->sc_r_power == 48)
1533 if (fep->sc_alpha > 0.1)
1535 gmx_fatal(FARGS, "sc_alpha (%f) for sc_r_power = 48 should usually be between 0.001 and 0.004", fep->sc_alpha);
1539 expand = ir->expandedvals;
1540 /* now read in the weights */
1541 parse_n_real(weights, &nweights, &(expand->init_lambda_weights));
1544 snew(expand->init_lambda_weights, fep->n_lambda); /* initialize to zero */
1546 else if (nweights != fep->n_lambda)
1548 gmx_fatal(FARGS, "Number of weights (%d) is not equal to number of lambda values (%d)",
1549 nweights, fep->n_lambda);
1551 if ((expand->nstexpanded < 0) && (ir->efep != efepNO))
1553 expand->nstexpanded = fep->nstdhdl;
1554 /* if you don't specify nstexpanded when doing expanded ensemble free energy calcs, it is set to nstdhdl */
1556 if ((expand->nstexpanded < 0) && ir->bSimTemp)
1558 expand->nstexpanded = 2*(int)(ir->opts.tau_t[0]/ir->delta_t);
1559 /* if you don't specify nstexpanded when doing expanded ensemble simulated tempering, it is set to
1560 2*tau_t just to be careful so it's not to frequent */
1565 static void do_simtemp_params(t_inputrec *ir)
1568 snew(ir->simtempvals->temperatures, ir->fepvals->n_lambda);
1569 GetSimTemps(ir->fepvals->n_lambda, ir->simtempvals, ir->fepvals->all_lambda[efptTEMPERATURE]);
1574 static void do_wall_params(t_inputrec *ir,
1575 char *wall_atomtype, char *wall_density,
1579 char *names[MAXPTR];
1582 opts->wall_atomtype[0] = NULL;
1583 opts->wall_atomtype[1] = NULL;
1585 ir->wall_atomtype[0] = -1;
1586 ir->wall_atomtype[1] = -1;
1587 ir->wall_density[0] = 0;
1588 ir->wall_density[1] = 0;
1592 nstr = str_nelem(wall_atomtype, MAXPTR, names);
1593 if (nstr != ir->nwall)
1595 gmx_fatal(FARGS, "Expected %d elements for wall_atomtype, found %d",
1598 for (i = 0; i < ir->nwall; i++)
1600 opts->wall_atomtype[i] = strdup(names[i]);
1603 if (ir->wall_type == ewt93 || ir->wall_type == ewt104)
1605 nstr = str_nelem(wall_density, MAXPTR, names);
1606 if (nstr != ir->nwall)
1608 gmx_fatal(FARGS, "Expected %d elements for wall-density, found %d", ir->nwall, nstr);
1610 for (i = 0; i < ir->nwall; i++)
1612 sscanf(names[i], "%lf", &dbl);
1615 gmx_fatal(FARGS, "wall-density[%d] = %f\n", i, dbl);
1617 ir->wall_density[i] = dbl;
1623 static void add_wall_energrps(gmx_groups_t *groups, int nwall, t_symtab *symtab)
1631 srenew(groups->grpname, groups->ngrpname+nwall);
1632 grps = &(groups->grps[egcENER]);
1633 srenew(grps->nm_ind, grps->nr+nwall);
1634 for (i = 0; i < nwall; i++)
1636 sprintf(str, "wall%d", i);
1637 groups->grpname[groups->ngrpname] = put_symtab(symtab, str);
1638 grps->nm_ind[grps->nr++] = groups->ngrpname++;
1643 void read_expandedparams(int *ninp_p, t_inpfile **inp_p,
1644 t_expanded *expand, warninp_t wi)
1646 int ninp, nerror = 0;
1652 /* read expanded ensemble parameters */
1653 CCTYPE ("expanded ensemble variables");
1654 ITYPE ("nstexpanded", expand->nstexpanded, -1);
1655 EETYPE("lmc-stats", expand->elamstats, elamstats_names);
1656 EETYPE("lmc-move", expand->elmcmove, elmcmove_names);
1657 EETYPE("lmc-weights-equil", expand->elmceq, elmceq_names);
1658 ITYPE ("weight-equil-number-all-lambda", expand->equil_n_at_lam, -1);
1659 ITYPE ("weight-equil-number-samples", expand->equil_samples, -1);
1660 ITYPE ("weight-equil-number-steps", expand->equil_steps, -1);
1661 RTYPE ("weight-equil-wl-delta", expand->equil_wl_delta, -1);
1662 RTYPE ("weight-equil-count-ratio", expand->equil_ratio, -1);
1663 CCTYPE("Seed for Monte Carlo in lambda space");
1664 ITYPE ("lmc-seed", expand->lmc_seed, -1);
1665 RTYPE ("mc-temperature", expand->mc_temp, -1);
1666 ITYPE ("lmc-repeats", expand->lmc_repeats, 1);
1667 ITYPE ("lmc-gibbsdelta", expand->gibbsdeltalam, -1);
1668 ITYPE ("lmc-forced-nstart", expand->lmc_forced_nstart, 0);
1669 EETYPE("symmetrized-transition-matrix", expand->bSymmetrizedTMatrix, yesno_names);
1670 ITYPE("nst-transition-matrix", expand->nstTij, -1);
1671 ITYPE ("mininum-var-min", expand->minvarmin, 100); /*default is reasonable */
1672 ITYPE ("weight-c-range", expand->c_range, 0); /* default is just C=0 */
1673 RTYPE ("wl-scale", expand->wl_scale, 0.8);
1674 RTYPE ("wl-ratio", expand->wl_ratio, 0.8);
1675 RTYPE ("init-wl-delta", expand->init_wl_delta, 1.0);
1676 EETYPE("wl-oneovert", expand->bWLoneovert, yesno_names);
1684 void get_ir(const char *mdparin, const char *mdparout,
1685 t_inputrec *ir, t_gromppopts *opts,
1689 double dumdub[2][6];
1693 char warn_buf[STRLEN];
1694 t_lambda *fep = ir->fepvals;
1695 t_expanded *expand = ir->expandedvals;
1697 inp = read_inpfile(mdparin, &ninp, wi);
1699 snew(dumstr[0], STRLEN);
1700 snew(dumstr[1], STRLEN);
1702 /* remove the following deprecated commands */
1705 REM_TYPE("domain-decomposition");
1706 REM_TYPE("andersen-seed");
1708 REM_TYPE("dihre-fc");
1709 REM_TYPE("dihre-tau");
1710 REM_TYPE("nstdihreout");
1711 REM_TYPE("nstcheckpoint");
1713 /* replace the following commands with the clearer new versions*/
1714 REPL_TYPE("unconstrained-start", "continuation");
1715 REPL_TYPE("foreign-lambda", "fep-lambdas");
1716 REPL_TYPE("verlet-buffer-drift", "verlet-buffer-tolerance");
1718 CCTYPE ("VARIOUS PREPROCESSING OPTIONS");
1719 CTYPE ("Preprocessor information: use cpp syntax.");
1720 CTYPE ("e.g.: -I/home/joe/doe -I/home/mary/roe");
1721 STYPE ("include", opts->include, NULL);
1722 CTYPE ("e.g.: -DPOSRES -DFLEXIBLE (note these variable names are case sensitive)");
1723 STYPE ("define", opts->define, NULL);
1725 CCTYPE ("RUN CONTROL PARAMETERS");
1726 EETYPE("integrator", ir->eI, ei_names);
1727 CTYPE ("Start time and timestep in ps");
1728 RTYPE ("tinit", ir->init_t, 0.0);
1729 RTYPE ("dt", ir->delta_t, 0.001);
1730 STEPTYPE ("nsteps", ir->nsteps, 0);
1731 CTYPE ("For exact run continuation or redoing part of a run");
1732 STEPTYPE ("init-step", ir->init_step, 0);
1733 CTYPE ("Part index is updated automatically on checkpointing (keeps files separate)");
1734 ITYPE ("simulation-part", ir->simulation_part, 1);
1735 CTYPE ("mode for center of mass motion removal");
1736 EETYPE("comm-mode", ir->comm_mode, ecm_names);
1737 CTYPE ("number of steps for center of mass motion removal");
1738 ITYPE ("nstcomm", ir->nstcomm, 100);
1739 CTYPE ("group(s) for center of mass motion removal");
1740 STYPE ("comm-grps", vcm, NULL);
1742 CCTYPE ("LANGEVIN DYNAMICS OPTIONS");
1743 CTYPE ("Friction coefficient (amu/ps) and random seed");
1744 RTYPE ("bd-fric", ir->bd_fric, 0.0);
1745 ITYPE ("ld-seed", ir->ld_seed, 1993);
1748 CCTYPE ("ENERGY MINIMIZATION OPTIONS");
1749 CTYPE ("Force tolerance and initial step-size");
1750 RTYPE ("emtol", ir->em_tol, 10.0);
1751 RTYPE ("emstep", ir->em_stepsize, 0.01);
1752 CTYPE ("Max number of iterations in relax-shells");
1753 ITYPE ("niter", ir->niter, 20);
1754 CTYPE ("Step size (ps^2) for minimization of flexible constraints");
1755 RTYPE ("fcstep", ir->fc_stepsize, 0);
1756 CTYPE ("Frequency of steepest descents steps when doing CG");
1757 ITYPE ("nstcgsteep", ir->nstcgsteep, 1000);
1758 ITYPE ("nbfgscorr", ir->nbfgscorr, 10);
1760 CCTYPE ("TEST PARTICLE INSERTION OPTIONS");
1761 RTYPE ("rtpi", ir->rtpi, 0.05);
1763 /* Output options */
1764 CCTYPE ("OUTPUT CONTROL OPTIONS");
1765 CTYPE ("Output frequency for coords (x), velocities (v) and forces (f)");
1766 ITYPE ("nstxout", ir->nstxout, 0);
1767 ITYPE ("nstvout", ir->nstvout, 0);
1768 ITYPE ("nstfout", ir->nstfout, 0);
1769 ir->nstcheckpoint = 1000;
1770 CTYPE ("Output frequency for energies to log file and energy file");
1771 ITYPE ("nstlog", ir->nstlog, 1000);
1772 ITYPE ("nstcalcenergy", ir->nstcalcenergy, 100);
1773 ITYPE ("nstenergy", ir->nstenergy, 1000);
1774 CTYPE ("Output frequency and precision for .xtc file");
1775 ITYPE ("nstxtcout", ir->nstxtcout, 0);
1776 RTYPE ("xtc-precision", ir->xtcprec, 1000.0);
1777 CTYPE ("This selects the subset of atoms for the .xtc file. You can");
1778 CTYPE ("select multiple groups. By default all atoms will be written.");
1779 STYPE ("xtc-grps", xtc_grps, NULL);
1780 CTYPE ("Selection of energy groups");
1781 STYPE ("energygrps", energy, NULL);
1783 /* Neighbor searching */
1784 CCTYPE ("NEIGHBORSEARCHING PARAMETERS");
1785 CTYPE ("cut-off scheme (Verlet: particle based cut-offs, group: using charge groups)");
1786 EETYPE("cutoff-scheme", ir->cutoff_scheme, ecutscheme_names);
1787 CTYPE ("nblist update frequency");
1788 ITYPE ("nstlist", ir->nstlist, 10);
1789 CTYPE ("ns algorithm (simple or grid)");
1790 EETYPE("ns-type", ir->ns_type, ens_names);
1791 /* set ndelta to the optimal value of 2 */
1793 CTYPE ("Periodic boundary conditions: xyz, no, xy");
1794 EETYPE("pbc", ir->ePBC, epbc_names);
1795 EETYPE("periodic-molecules", ir->bPeriodicMols, yesno_names);
1796 CTYPE ("Allowed energy error due to the Verlet buffer in kJ/mol/ps per atom,");
1797 CTYPE ("a value of -1 means: use rlist");
1798 RTYPE("verlet-buffer-tolerance", ir->verletbuf_tol, 0.005);
1799 CTYPE ("nblist cut-off");
1800 RTYPE ("rlist", ir->rlist, 1.0);
1801 CTYPE ("long-range cut-off for switched potentials");
1802 RTYPE ("rlistlong", ir->rlistlong, -1);
1803 ITYPE ("nstcalclr", ir->nstcalclr, -1);
1805 /* Electrostatics */
1806 CCTYPE ("OPTIONS FOR ELECTROSTATICS AND VDW");
1807 CTYPE ("Method for doing electrostatics");
1808 EETYPE("coulombtype", ir->coulombtype, eel_names);
1809 EETYPE("coulomb-modifier", ir->coulomb_modifier, eintmod_names);
1810 CTYPE ("cut-off lengths");
1811 RTYPE ("rcoulomb-switch", ir->rcoulomb_switch, 0.0);
1812 RTYPE ("rcoulomb", ir->rcoulomb, 1.0);
1813 CTYPE ("Relative dielectric constant for the medium and the reaction field");
1814 RTYPE ("epsilon-r", ir->epsilon_r, 1.0);
1815 RTYPE ("epsilon-rf", ir->epsilon_rf, 0.0);
1816 CTYPE ("Method for doing Van der Waals");
1817 EETYPE("vdw-type", ir->vdwtype, evdw_names);
1818 EETYPE("vdw-modifier", ir->vdw_modifier, eintmod_names);
1819 CTYPE ("cut-off lengths");
1820 RTYPE ("rvdw-switch", ir->rvdw_switch, 0.0);
1821 RTYPE ("rvdw", ir->rvdw, 1.0);
1822 CTYPE ("Apply long range dispersion corrections for Energy and Pressure");
1823 EETYPE("DispCorr", ir->eDispCorr, edispc_names);
1824 CTYPE ("Extension of the potential lookup tables beyond the cut-off");
1825 RTYPE ("table-extension", ir->tabext, 1.0);
1826 CTYPE ("Separate tables between energy group pairs");
1827 STYPE ("energygrp-table", egptable, NULL);
1828 CTYPE ("Spacing for the PME/PPPM FFT grid");
1829 RTYPE ("fourierspacing", ir->fourier_spacing, 0.12);
1830 CTYPE ("FFT grid size, when a value is 0 fourierspacing will be used");
1831 ITYPE ("fourier-nx", ir->nkx, 0);
1832 ITYPE ("fourier-ny", ir->nky, 0);
1833 ITYPE ("fourier-nz", ir->nkz, 0);
1834 CTYPE ("EWALD/PME/PPPM parameters");
1835 ITYPE ("pme-order", ir->pme_order, 4);
1836 RTYPE ("ewald-rtol", ir->ewald_rtol, 0.00001);
1837 RTYPE ("ewald-rtol-lj", ir->ewald_rtol_lj, 0.001);
1838 EETYPE("lj-pme-comb-rule", ir->ljpme_combination_rule, eljpme_names);
1839 EETYPE("ewald-geometry", ir->ewald_geometry, eewg_names);
1840 RTYPE ("epsilon-surface", ir->epsilon_surface, 0.0);
1841 EETYPE("optimize-fft", ir->bOptFFT, yesno_names);
1843 CCTYPE("IMPLICIT SOLVENT ALGORITHM");
1844 EETYPE("implicit-solvent", ir->implicit_solvent, eis_names);
1846 CCTYPE ("GENERALIZED BORN ELECTROSTATICS");
1847 CTYPE ("Algorithm for calculating Born radii");
1848 EETYPE("gb-algorithm", ir->gb_algorithm, egb_names);
1849 CTYPE ("Frequency of calculating the Born radii inside rlist");
1850 ITYPE ("nstgbradii", ir->nstgbradii, 1);
1851 CTYPE ("Cutoff for Born radii calculation; the contribution from atoms");
1852 CTYPE ("between rlist and rgbradii is updated every nstlist steps");
1853 RTYPE ("rgbradii", ir->rgbradii, 1.0);
1854 CTYPE ("Dielectric coefficient of the implicit solvent");
1855 RTYPE ("gb-epsilon-solvent", ir->gb_epsilon_solvent, 80.0);
1856 CTYPE ("Salt concentration in M for Generalized Born models");
1857 RTYPE ("gb-saltconc", ir->gb_saltconc, 0.0);
1858 CTYPE ("Scaling factors used in the OBC GB model. Default values are OBC(II)");
1859 RTYPE ("gb-obc-alpha", ir->gb_obc_alpha, 1.0);
1860 RTYPE ("gb-obc-beta", ir->gb_obc_beta, 0.8);
1861 RTYPE ("gb-obc-gamma", ir->gb_obc_gamma, 4.85);
1862 RTYPE ("gb-dielectric-offset", ir->gb_dielectric_offset, 0.009);
1863 EETYPE("sa-algorithm", ir->sa_algorithm, esa_names);
1864 CTYPE ("Surface tension (kJ/mol/nm^2) for the SA (nonpolar surface) part of GBSA");
1865 CTYPE ("The value -1 will set default value for Still/HCT/OBC GB-models.");
1866 RTYPE ("sa-surface-tension", ir->sa_surface_tension, -1);
1868 /* Coupling stuff */
1869 CCTYPE ("OPTIONS FOR WEAK COUPLING ALGORITHMS");
1870 CTYPE ("Temperature coupling");
1871 EETYPE("tcoupl", ir->etc, etcoupl_names);
1872 ITYPE ("nsttcouple", ir->nsttcouple, -1);
1873 ITYPE("nh-chain-length", ir->opts.nhchainlength, NHCHAINLENGTH);
1874 EETYPE("print-nose-hoover-chain-variables", ir->bPrintNHChains, yesno_names);
1875 CTYPE ("Groups to couple separately");
1876 STYPE ("tc-grps", tcgrps, NULL);
1877 CTYPE ("Time constant (ps) and reference temperature (K)");
1878 STYPE ("tau-t", tau_t, NULL);
1879 STYPE ("ref-t", ref_t, NULL);
1880 CTYPE ("pressure coupling");
1881 EETYPE("pcoupl", ir->epc, epcoupl_names);
1882 EETYPE("pcoupltype", ir->epct, epcoupltype_names);
1883 ITYPE ("nstpcouple", ir->nstpcouple, -1);
1884 CTYPE ("Time constant (ps), compressibility (1/bar) and reference P (bar)");
1885 RTYPE ("tau-p", ir->tau_p, 1.0);
1886 STYPE ("compressibility", dumstr[0], NULL);
1887 STYPE ("ref-p", dumstr[1], NULL);
1888 CTYPE ("Scaling of reference coordinates, No, All or COM");
1889 EETYPE ("refcoord-scaling", ir->refcoord_scaling, erefscaling_names);
1892 CCTYPE ("OPTIONS FOR QMMM calculations");
1893 EETYPE("QMMM", ir->bQMMM, yesno_names);
1894 CTYPE ("Groups treated Quantum Mechanically");
1895 STYPE ("QMMM-grps", QMMM, NULL);
1896 CTYPE ("QM method");
1897 STYPE("QMmethod", QMmethod, NULL);
1898 CTYPE ("QMMM scheme");
1899 EETYPE("QMMMscheme", ir->QMMMscheme, eQMMMscheme_names);
1900 CTYPE ("QM basisset");
1901 STYPE("QMbasis", QMbasis, NULL);
1902 CTYPE ("QM charge");
1903 STYPE ("QMcharge", QMcharge, NULL);
1904 CTYPE ("QM multiplicity");
1905 STYPE ("QMmult", QMmult, NULL);
1906 CTYPE ("Surface Hopping");
1907 STYPE ("SH", bSH, NULL);
1908 CTYPE ("CAS space options");
1909 STYPE ("CASorbitals", CASorbitals, NULL);
1910 STYPE ("CASelectrons", CASelectrons, NULL);
1911 STYPE ("SAon", SAon, NULL);
1912 STYPE ("SAoff", SAoff, NULL);
1913 STYPE ("SAsteps", SAsteps, NULL);
1914 CTYPE ("Scale factor for MM charges");
1915 RTYPE ("MMChargeScaleFactor", ir->scalefactor, 1.0);
1916 CTYPE ("Optimization of QM subsystem");
1917 STYPE ("bOPT", bOPT, NULL);
1918 STYPE ("bTS", bTS, NULL);
1920 /* Simulated annealing */
1921 CCTYPE("SIMULATED ANNEALING");
1922 CTYPE ("Type of annealing for each temperature group (no/single/periodic)");
1923 STYPE ("annealing", anneal, NULL);
1924 CTYPE ("Number of time points to use for specifying annealing in each group");
1925 STYPE ("annealing-npoints", anneal_npoints, NULL);
1926 CTYPE ("List of times at the annealing points for each group");
1927 STYPE ("annealing-time", anneal_time, NULL);
1928 CTYPE ("Temp. at each annealing point, for each group.");
1929 STYPE ("annealing-temp", anneal_temp, NULL);
1932 CCTYPE ("GENERATE VELOCITIES FOR STARTUP RUN");
1933 EETYPE("gen-vel", opts->bGenVel, yesno_names);
1934 RTYPE ("gen-temp", opts->tempi, 300.0);
1935 ITYPE ("gen-seed", opts->seed, 173529);
1938 CCTYPE ("OPTIONS FOR BONDS");
1939 EETYPE("constraints", opts->nshake, constraints);
1940 CTYPE ("Type of constraint algorithm");
1941 EETYPE("constraint-algorithm", ir->eConstrAlg, econstr_names);
1942 CTYPE ("Do not constrain the start configuration");
1943 EETYPE("continuation", ir->bContinuation, yesno_names);
1944 CTYPE ("Use successive overrelaxation to reduce the number of shake iterations");
1945 EETYPE("Shake-SOR", ir->bShakeSOR, yesno_names);
1946 CTYPE ("Relative tolerance of shake");
1947 RTYPE ("shake-tol", ir->shake_tol, 0.0001);
1948 CTYPE ("Highest order in the expansion of the constraint coupling matrix");
1949 ITYPE ("lincs-order", ir->nProjOrder, 4);
1950 CTYPE ("Number of iterations in the final step of LINCS. 1 is fine for");
1951 CTYPE ("normal simulations, but use 2 to conserve energy in NVE runs.");
1952 CTYPE ("For energy minimization with constraints it should be 4 to 8.");
1953 ITYPE ("lincs-iter", ir->nLincsIter, 1);
1954 CTYPE ("Lincs will write a warning to the stderr if in one step a bond");
1955 CTYPE ("rotates over more degrees than");
1956 RTYPE ("lincs-warnangle", ir->LincsWarnAngle, 30.0);
1957 CTYPE ("Convert harmonic bonds to morse potentials");
1958 EETYPE("morse", opts->bMorse, yesno_names);
1960 /* Energy group exclusions */
1961 CCTYPE ("ENERGY GROUP EXCLUSIONS");
1962 CTYPE ("Pairs of energy groups for which all non-bonded interactions are excluded");
1963 STYPE ("energygrp-excl", egpexcl, NULL);
1967 CTYPE ("Number of walls, type, atom types, densities and box-z scale factor for Ewald");
1968 ITYPE ("nwall", ir->nwall, 0);
1969 EETYPE("wall-type", ir->wall_type, ewt_names);
1970 RTYPE ("wall-r-linpot", ir->wall_r_linpot, -1);
1971 STYPE ("wall-atomtype", wall_atomtype, NULL);
1972 STYPE ("wall-density", wall_density, NULL);
1973 RTYPE ("wall-ewald-zfac", ir->wall_ewald_zfac, 3);
1976 CCTYPE("COM PULLING");
1977 CTYPE("Pull type: no, umbrella, constraint or constant-force");
1978 EETYPE("pull", ir->ePull, epull_names);
1979 if (ir->ePull != epullNO)
1982 pull_grp = read_pullparams(&ninp, &inp, ir->pull, &opts->pull_start, wi);
1985 /* Enforced rotation */
1986 CCTYPE("ENFORCED ROTATION");
1987 CTYPE("Enforced rotation: No or Yes");
1988 EETYPE("rotation", ir->bRot, yesno_names);
1992 rot_grp = read_rotparams(&ninp, &inp, ir->rot, wi);
1996 CCTYPE("NMR refinement stuff");
1997 CTYPE ("Distance restraints type: No, Simple or Ensemble");
1998 EETYPE("disre", ir->eDisre, edisre_names);
1999 CTYPE ("Force weighting of pairs in one distance restraint: Conservative or Equal");
2000 EETYPE("disre-weighting", ir->eDisreWeighting, edisreweighting_names);
2001 CTYPE ("Use sqrt of the time averaged times the instantaneous violation");
2002 EETYPE("disre-mixed", ir->bDisreMixed, yesno_names);
2003 RTYPE ("disre-fc", ir->dr_fc, 1000.0);
2004 RTYPE ("disre-tau", ir->dr_tau, 0.0);
2005 CTYPE ("Output frequency for pair distances to energy file");
2006 ITYPE ("nstdisreout", ir->nstdisreout, 100);
2007 CTYPE ("Orientation restraints: No or Yes");
2008 EETYPE("orire", opts->bOrire, yesno_names);
2009 CTYPE ("Orientation restraints force constant and tau for time averaging");
2010 RTYPE ("orire-fc", ir->orires_fc, 0.0);
2011 RTYPE ("orire-tau", ir->orires_tau, 0.0);
2012 STYPE ("orire-fitgrp", orirefitgrp, NULL);
2013 CTYPE ("Output frequency for trace(SD) and S to energy file");
2014 ITYPE ("nstorireout", ir->nstorireout, 100);
2016 /* free energy variables */
2017 CCTYPE ("Free energy variables");
2018 EETYPE("free-energy", ir->efep, efep_names);
2019 STYPE ("couple-moltype", couple_moltype, NULL);
2020 EETYPE("couple-lambda0", opts->couple_lam0, couple_lam);
2021 EETYPE("couple-lambda1", opts->couple_lam1, couple_lam);
2022 EETYPE("couple-intramol", opts->bCoupleIntra, yesno_names);
2024 RTYPE ("init-lambda", fep->init_lambda, -1); /* start with -1 so
2026 it was not entered */
2027 ITYPE ("init-lambda-state", fep->init_fep_state, -1);
2028 RTYPE ("delta-lambda", fep->delta_lambda, 0.0);
2029 ITYPE ("nstdhdl", fep->nstdhdl, 50);
2030 STYPE ("fep-lambdas", fep_lambda[efptFEP], NULL);
2031 STYPE ("mass-lambdas", fep_lambda[efptMASS], NULL);
2032 STYPE ("coul-lambdas", fep_lambda[efptCOUL], NULL);
2033 STYPE ("vdw-lambdas", fep_lambda[efptVDW], NULL);
2034 STYPE ("bonded-lambdas", fep_lambda[efptBONDED], NULL);
2035 STYPE ("restraint-lambdas", fep_lambda[efptRESTRAINT], NULL);
2036 STYPE ("temperature-lambdas", fep_lambda[efptTEMPERATURE], NULL);
2037 ITYPE ("calc-lambda-neighbors", fep->lambda_neighbors, 1);
2038 STYPE ("init-lambda-weights", lambda_weights, NULL);
2039 EETYPE("dhdl-print-energy", fep->bPrintEnergy, yesno_names);
2040 RTYPE ("sc-alpha", fep->sc_alpha, 0.0);
2041 ITYPE ("sc-power", fep->sc_power, 1);
2042 RTYPE ("sc-r-power", fep->sc_r_power, 6.0);
2043 RTYPE ("sc-sigma", fep->sc_sigma, 0.3);
2044 EETYPE("sc-coul", fep->bScCoul, yesno_names);
2045 ITYPE ("dh_hist_size", fep->dh_hist_size, 0);
2046 RTYPE ("dh_hist_spacing", fep->dh_hist_spacing, 0.1);
2047 EETYPE("separate-dhdl-file", fep->separate_dhdl_file,
2048 separate_dhdl_file_names);
2049 EETYPE("dhdl-derivatives", fep->dhdl_derivatives, dhdl_derivatives_names);
2050 ITYPE ("dh_hist_size", fep->dh_hist_size, 0);
2051 RTYPE ("dh_hist_spacing", fep->dh_hist_spacing, 0.1);
2053 /* Non-equilibrium MD stuff */
2054 CCTYPE("Non-equilibrium MD stuff");
2055 STYPE ("acc-grps", accgrps, NULL);
2056 STYPE ("accelerate", acc, NULL);
2057 STYPE ("freezegrps", freeze, NULL);
2058 STYPE ("freezedim", frdim, NULL);
2059 RTYPE ("cos-acceleration", ir->cos_accel, 0);
2060 STYPE ("deform", deform, NULL);
2062 /* simulated tempering variables */
2063 CCTYPE("simulated tempering variables");
2064 EETYPE("simulated-tempering", ir->bSimTemp, yesno_names);
2065 EETYPE("simulated-tempering-scaling", ir->simtempvals->eSimTempScale, esimtemp_names);
2066 RTYPE("sim-temp-low", ir->simtempvals->simtemp_low, 300.0);
2067 RTYPE("sim-temp-high", ir->simtempvals->simtemp_high, 300.0);
2069 /* expanded ensemble variables */
2070 if (ir->efep == efepEXPANDED || ir->bSimTemp)
2072 read_expandedparams(&ninp, &inp, expand, wi);
2075 /* Electric fields */
2076 CCTYPE("Electric fields");
2077 CTYPE ("Format is number of terms (int) and for all terms an amplitude (real)");
2078 CTYPE ("and a phase angle (real)");
2079 STYPE ("E-x", efield_x, NULL);
2080 STYPE ("E-xt", efield_xt, NULL);
2081 STYPE ("E-y", efield_y, NULL);
2082 STYPE ("E-yt", efield_yt, NULL);
2083 STYPE ("E-z", efield_z, NULL);
2084 STYPE ("E-zt", efield_zt, NULL);
2086 /* AdResS defined thingies */
2087 CCTYPE ("AdResS parameters");
2088 EETYPE("adress", ir->bAdress, yesno_names);
2091 snew(ir->adress, 1);
2092 read_adressparams(&ninp, &inp, ir->adress, wi);
2095 /* User defined thingies */
2096 CCTYPE ("User defined thingies");
2097 STYPE ("user1-grps", user1, NULL);
2098 STYPE ("user2-grps", user2, NULL);
2099 ITYPE ("userint1", ir->userint1, 0);
2100 ITYPE ("userint2", ir->userint2, 0);
2101 ITYPE ("userint3", ir->userint3, 0);
2102 ITYPE ("userint4", ir->userint4, 0);
2103 RTYPE ("userreal1", ir->userreal1, 0);
2104 RTYPE ("userreal2", ir->userreal2, 0);
2105 RTYPE ("userreal3", ir->userreal3, 0);
2106 RTYPE ("userreal4", ir->userreal4, 0);
2109 write_inpfile(mdparout, ninp, inp, FALSE, wi);
2110 for (i = 0; (i < ninp); i++)
2113 sfree(inp[i].value);
2117 /* Process options if necessary */
2118 for (m = 0; m < 2; m++)
2120 for (i = 0; i < 2*DIM; i++)
2129 if (sscanf(dumstr[m], "%lf", &(dumdub[m][XX])) != 1)
2131 warning_error(wi, "Pressure coupling not enough values (I need 1)");
2133 dumdub[m][YY] = dumdub[m][ZZ] = dumdub[m][XX];
2135 case epctSEMIISOTROPIC:
2136 case epctSURFACETENSION:
2137 if (sscanf(dumstr[m], "%lf%lf",
2138 &(dumdub[m][XX]), &(dumdub[m][ZZ])) != 2)
2140 warning_error(wi, "Pressure coupling not enough values (I need 2)");
2142 dumdub[m][YY] = dumdub[m][XX];
2144 case epctANISOTROPIC:
2145 if (sscanf(dumstr[m], "%lf%lf%lf%lf%lf%lf",
2146 &(dumdub[m][XX]), &(dumdub[m][YY]), &(dumdub[m][ZZ]),
2147 &(dumdub[m][3]), &(dumdub[m][4]), &(dumdub[m][5])) != 6)
2149 warning_error(wi, "Pressure coupling not enough values (I need 6)");
2153 gmx_fatal(FARGS, "Pressure coupling type %s not implemented yet",
2154 epcoupltype_names[ir->epct]);
2158 clear_mat(ir->ref_p);
2159 clear_mat(ir->compress);
2160 for (i = 0; i < DIM; i++)
2162 ir->ref_p[i][i] = dumdub[1][i];
2163 ir->compress[i][i] = dumdub[0][i];
2165 if (ir->epct == epctANISOTROPIC)
2167 ir->ref_p[XX][YY] = dumdub[1][3];
2168 ir->ref_p[XX][ZZ] = dumdub[1][4];
2169 ir->ref_p[YY][ZZ] = dumdub[1][5];
2170 if (ir->ref_p[XX][YY] != 0 && ir->ref_p[XX][ZZ] != 0 && ir->ref_p[YY][ZZ] != 0)
2172 warning(wi, "All off-diagonal reference pressures are non-zero. Are you sure you want to apply a threefold shear stress?\n");
2174 ir->compress[XX][YY] = dumdub[0][3];
2175 ir->compress[XX][ZZ] = dumdub[0][4];
2176 ir->compress[YY][ZZ] = dumdub[0][5];
2177 for (i = 0; i < DIM; i++)
2179 for (m = 0; m < i; m++)
2181 ir->ref_p[i][m] = ir->ref_p[m][i];
2182 ir->compress[i][m] = ir->compress[m][i];
2187 if (ir->comm_mode == ecmNO)
2192 opts->couple_moltype = NULL;
2193 if (strlen(couple_moltype) > 0)
2195 if (ir->efep != efepNO)
2197 opts->couple_moltype = strdup(couple_moltype);
2198 if (opts->couple_lam0 == opts->couple_lam1)
2200 warning(wi, "The lambda=0 and lambda=1 states for coupling are identical");
2202 if (ir->eI == eiMD && (opts->couple_lam0 == ecouplamNONE ||
2203 opts->couple_lam1 == ecouplamNONE))
2205 warning(wi, "For proper sampling of the (nearly) decoupled state, stochastic dynamics should be used");
2210 warning(wi, "Can not couple a molecule with free_energy = no");
2213 /* FREE ENERGY AND EXPANDED ENSEMBLE OPTIONS */
2214 if (ir->efep != efepNO)
2216 if (fep->delta_lambda > 0)
2218 ir->efep = efepSLOWGROWTH;
2224 fep->bPrintEnergy = TRUE;
2225 /* always print out the energy to dhdl if we are doing expanded ensemble, since we need the total energy
2226 if the temperature is changing. */
2229 if ((ir->efep != efepNO) || ir->bSimTemp)
2231 ir->bExpanded = FALSE;
2232 if ((ir->efep == efepEXPANDED) || ir->bSimTemp)
2234 ir->bExpanded = TRUE;
2236 do_fep_params(ir, fep_lambda, lambda_weights);
2237 if (ir->bSimTemp) /* done after fep params */
2239 do_simtemp_params(ir);
2244 ir->fepvals->n_lambda = 0;
2247 /* WALL PARAMETERS */
2249 do_wall_params(ir, wall_atomtype, wall_density, opts);
2251 /* ORIENTATION RESTRAINT PARAMETERS */
2253 if (opts->bOrire && str_nelem(orirefitgrp, MAXPTR, NULL) != 1)
2255 warning_error(wi, "ERROR: Need one orientation restraint fit group\n");
2258 /* DEFORMATION PARAMETERS */
2260 clear_mat(ir->deform);
2261 for (i = 0; i < 6; i++)
2265 m = sscanf(deform, "%lf %lf %lf %lf %lf %lf",
2266 &(dumdub[0][0]), &(dumdub[0][1]), &(dumdub[0][2]),
2267 &(dumdub[0][3]), &(dumdub[0][4]), &(dumdub[0][5]));
2268 for (i = 0; i < 3; i++)
2270 ir->deform[i][i] = dumdub[0][i];
2272 ir->deform[YY][XX] = dumdub[0][3];
2273 ir->deform[ZZ][XX] = dumdub[0][4];
2274 ir->deform[ZZ][YY] = dumdub[0][5];
2275 if (ir->epc != epcNO)
2277 for (i = 0; i < 3; i++)
2279 for (j = 0; j <= i; j++)
2281 if (ir->deform[i][j] != 0 && ir->compress[i][j] != 0)
2283 warning_error(wi, "A box element has deform set and compressibility > 0");
2287 for (i = 0; i < 3; i++)
2289 for (j = 0; j < i; j++)
2291 if (ir->deform[i][j] != 0)
2293 for (m = j; m < DIM; m++)
2295 if (ir->compress[m][j] != 0)
2297 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.");
2298 warning(wi, warn_buf);
2310 static int search_QMstring(const char *s, int ng, const char *gn[])
2312 /* same as normal search_string, but this one searches QM strings */
2315 for (i = 0; (i < ng); i++)
2317 if (gmx_strcasecmp(s, gn[i]) == 0)
2323 gmx_fatal(FARGS, "this QM method or basisset (%s) is not implemented\n!", s);
2327 } /* search_QMstring */
2329 /* We would like gn to be const as well, but C doesn't allow this */
2330 int search_string(const char *s, int ng, char *gn[])
2334 for (i = 0; (i < ng); i++)
2336 if (gmx_strcasecmp(s, gn[i]) == 0)
2343 "Group %s referenced in the .mdp file was not found in the index file.\n"
2344 "Group names must match either [moleculetype] names or custom index group\n"
2345 "names, in which case you must supply an index file to the '-n' option\n"
2352 static gmx_bool do_numbering(int natoms, gmx_groups_t *groups, int ng, char *ptrs[],
2353 t_blocka *block, char *gnames[],
2354 int gtype, int restnm,
2355 int grptp, gmx_bool bVerbose,
2358 unsigned short *cbuf;
2359 t_grps *grps = &(groups->grps[gtype]);
2360 int i, j, gid, aj, ognr, ntot = 0;
2363 char warn_buf[STRLEN];
2367 fprintf(debug, "Starting numbering %d groups of type %d\n", ng, gtype);
2370 title = gtypes[gtype];
2373 /* Mark all id's as not set */
2374 for (i = 0; (i < natoms); i++)
2379 snew(grps->nm_ind, ng+1); /* +1 for possible rest group */
2380 for (i = 0; (i < ng); i++)
2382 /* Lookup the group name in the block structure */
2383 gid = search_string(ptrs[i], block->nr, gnames);
2384 if ((grptp != egrptpONE) || (i == 0))
2386 grps->nm_ind[grps->nr++] = gid;
2390 fprintf(debug, "Found gid %d for group %s\n", gid, ptrs[i]);
2393 /* Now go over the atoms in the group */
2394 for (j = block->index[gid]; (j < block->index[gid+1]); j++)
2399 /* Range checking */
2400 if ((aj < 0) || (aj >= natoms))
2402 gmx_fatal(FARGS, "Invalid atom number %d in indexfile", aj);
2404 /* Lookup up the old group number */
2408 gmx_fatal(FARGS, "Atom %d in multiple %s groups (%d and %d)",
2409 aj+1, title, ognr+1, i+1);
2413 /* Store the group number in buffer */
2414 if (grptp == egrptpONE)
2427 /* Now check whether we have done all atoms */
2431 if (grptp == egrptpALL)
2433 gmx_fatal(FARGS, "%d atoms are not part of any of the %s groups",
2434 natoms-ntot, title);
2436 else if (grptp == egrptpPART)
2438 sprintf(warn_buf, "%d atoms are not part of any of the %s groups",
2439 natoms-ntot, title);
2440 warning_note(wi, warn_buf);
2442 /* Assign all atoms currently unassigned to a rest group */
2443 for (j = 0; (j < natoms); j++)
2445 if (cbuf[j] == NOGID)
2451 if (grptp != egrptpPART)
2456 "Making dummy/rest group for %s containing %d elements\n",
2457 title, natoms-ntot);
2459 /* Add group name "rest" */
2460 grps->nm_ind[grps->nr] = restnm;
2462 /* Assign the rest name to all atoms not currently assigned to a group */
2463 for (j = 0; (j < natoms); j++)
2465 if (cbuf[j] == NOGID)
2474 if (grps->nr == 1 && (ntot == 0 || ntot == natoms))
2476 /* All atoms are part of one (or no) group, no index required */
2477 groups->ngrpnr[gtype] = 0;
2478 groups->grpnr[gtype] = NULL;
2482 groups->ngrpnr[gtype] = natoms;
2483 snew(groups->grpnr[gtype], natoms);
2484 for (j = 0; (j < natoms); j++)
2486 groups->grpnr[gtype][j] = cbuf[j];
2492 return (bRest && grptp == egrptpPART);
2495 static void calc_nrdf(gmx_mtop_t *mtop, t_inputrec *ir, char **gnames)
2498 gmx_groups_t *groups;
2500 int natoms, ai, aj, i, j, d, g, imin, jmin;
2502 int *nrdf2, *na_vcm, na_tot;
2503 double *nrdf_tc, *nrdf_vcm, nrdf_uc, n_sub = 0;
2504 gmx_mtop_atomloop_all_t aloop;
2506 int mb, mol, ftype, as;
2507 gmx_molblock_t *molb;
2508 gmx_moltype_t *molt;
2511 * First calc 3xnr-atoms for each group
2512 * then subtract half a degree of freedom for each constraint
2514 * Only atoms and nuclei contribute to the degrees of freedom...
2519 groups = &mtop->groups;
2520 natoms = mtop->natoms;
2522 /* Allocate one more for a possible rest group */
2523 /* We need to sum degrees of freedom into doubles,
2524 * since floats give too low nrdf's above 3 million atoms.
2526 snew(nrdf_tc, groups->grps[egcTC].nr+1);
2527 snew(nrdf_vcm, groups->grps[egcVCM].nr+1);
2528 snew(na_vcm, groups->grps[egcVCM].nr+1);
2530 for (i = 0; i < groups->grps[egcTC].nr; i++)
2534 for (i = 0; i < groups->grps[egcVCM].nr+1; i++)
2539 snew(nrdf2, natoms);
2540 aloop = gmx_mtop_atomloop_all_init(mtop);
2541 while (gmx_mtop_atomloop_all_next(aloop, &i, &atom))
2544 if (atom->ptype == eptAtom || atom->ptype == eptNucleus)
2546 g = ggrpnr(groups, egcFREEZE, i);
2547 /* Double count nrdf for particle i */
2548 for (d = 0; d < DIM; d++)
2550 if (opts->nFreeze[g][d] == 0)
2555 nrdf_tc [ggrpnr(groups, egcTC, i)] += 0.5*nrdf2[i];
2556 nrdf_vcm[ggrpnr(groups, egcVCM, i)] += 0.5*nrdf2[i];
2561 for (mb = 0; mb < mtop->nmolblock; mb++)
2563 molb = &mtop->molblock[mb];
2564 molt = &mtop->moltype[molb->type];
2565 atom = molt->atoms.atom;
2566 for (mol = 0; mol < molb->nmol; mol++)
2568 for (ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++)
2570 ia = molt->ilist[ftype].iatoms;
2571 for (i = 0; i < molt->ilist[ftype].nr; )
2573 /* Subtract degrees of freedom for the constraints,
2574 * if the particles still have degrees of freedom left.
2575 * If one of the particles is a vsite or a shell, then all
2576 * constraint motion will go there, but since they do not
2577 * contribute to the constraints the degrees of freedom do not
2582 if (((atom[ia[1]].ptype == eptNucleus) ||
2583 (atom[ia[1]].ptype == eptAtom)) &&
2584 ((atom[ia[2]].ptype == eptNucleus) ||
2585 (atom[ia[2]].ptype == eptAtom)))
2603 imin = min(imin, nrdf2[ai]);
2604 jmin = min(jmin, nrdf2[aj]);
2607 nrdf_tc [ggrpnr(groups, egcTC, ai)] -= 0.5*imin;
2608 nrdf_tc [ggrpnr(groups, egcTC, aj)] -= 0.5*jmin;
2609 nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5*imin;
2610 nrdf_vcm[ggrpnr(groups, egcVCM, aj)] -= 0.5*jmin;
2612 ia += interaction_function[ftype].nratoms+1;
2613 i += interaction_function[ftype].nratoms+1;
2616 ia = molt->ilist[F_SETTLE].iatoms;
2617 for (i = 0; i < molt->ilist[F_SETTLE].nr; )
2619 /* Subtract 1 dof from every atom in the SETTLE */
2620 for (j = 0; j < 3; j++)
2623 imin = min(2, nrdf2[ai]);
2625 nrdf_tc [ggrpnr(groups, egcTC, ai)] -= 0.5*imin;
2626 nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5*imin;
2631 as += molt->atoms.nr;
2635 if (ir->ePull == epullCONSTRAINT)
2637 /* Correct nrdf for the COM constraints.
2638 * We correct using the TC and VCM group of the first atom
2639 * in the reference and pull group. If atoms in one pull group
2640 * belong to different TC or VCM groups it is anyhow difficult
2641 * to determine the optimal nrdf assignment.
2645 for (i = 0; i < pull->ncoord; i++)
2649 for (j = 0; j < 2; j++)
2651 const t_pull_group *pgrp;
2653 pgrp = &pull->group[pull->coord[i].group[j]];
2657 /* Subtract 1/2 dof from each group */
2659 nrdf_tc [ggrpnr(groups, egcTC, ai)] -= 0.5*imin;
2660 nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5*imin;
2661 if (nrdf_tc[ggrpnr(groups, egcTC, ai)] < 0)
2663 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)]]);
2668 /* We need to subtract the whole DOF from group j=1 */
2675 if (ir->nstcomm != 0)
2677 /* Subtract 3 from the number of degrees of freedom in each vcm group
2678 * when com translation is removed and 6 when rotation is removed
2681 switch (ir->comm_mode)
2684 n_sub = ndof_com(ir);
2691 gmx_incons("Checking comm_mode");
2694 for (i = 0; i < groups->grps[egcTC].nr; i++)
2696 /* Count the number of atoms of TC group i for every VCM group */
2697 for (j = 0; j < groups->grps[egcVCM].nr+1; j++)
2702 for (ai = 0; ai < natoms; ai++)
2704 if (ggrpnr(groups, egcTC, ai) == i)
2706 na_vcm[ggrpnr(groups, egcVCM, ai)]++;
2710 /* Correct for VCM removal according to the fraction of each VCM
2711 * group present in this TC group.
2713 nrdf_uc = nrdf_tc[i];
2716 fprintf(debug, "T-group[%d] nrdf_uc = %g, n_sub = %g\n",
2720 for (j = 0; j < groups->grps[egcVCM].nr+1; j++)
2722 if (nrdf_vcm[j] > n_sub)
2724 nrdf_tc[i] += nrdf_uc*((double)na_vcm[j]/(double)na_tot)*
2725 (nrdf_vcm[j] - n_sub)/nrdf_vcm[j];
2729 fprintf(debug, " nrdf_vcm[%d] = %g, nrdf = %g\n",
2730 j, nrdf_vcm[j], nrdf_tc[i]);
2735 for (i = 0; (i < groups->grps[egcTC].nr); i++)
2737 opts->nrdf[i] = nrdf_tc[i];
2738 if (opts->nrdf[i] < 0)
2743 "Number of degrees of freedom in T-Coupling group %s is %.2f\n",
2744 gnames[groups->grps[egcTC].nm_ind[i]], opts->nrdf[i]);
2753 static void decode_cos(char *s, t_cosines *cosine)
2756 char format[STRLEN], f1[STRLEN];
2768 sscanf(t, "%d", &(cosine->n));
2775 snew(cosine->a, cosine->n);
2776 snew(cosine->phi, cosine->n);
2778 sprintf(format, "%%*d");
2779 for (i = 0; (i < cosine->n); i++)
2782 strcat(f1, "%lf%lf");
2783 if (sscanf(t, f1, &a, &phi) < 2)
2785 gmx_fatal(FARGS, "Invalid input for electric field shift: '%s'", t);
2788 cosine->phi[i] = phi;
2789 strcat(format, "%*lf%*lf");
2796 static gmx_bool do_egp_flag(t_inputrec *ir, gmx_groups_t *groups,
2797 const char *option, const char *val, int flag)
2799 /* The maximum number of energy group pairs would be MAXPTR*(MAXPTR+1)/2.
2800 * But since this is much larger than STRLEN, such a line can not be parsed.
2801 * The real maximum is the number of names that fit in a string: STRLEN/2.
2803 #define EGP_MAX (STRLEN/2)
2804 int nelem, i, j, k, nr;
2805 char *names[EGP_MAX];
2809 gnames = groups->grpname;
2811 nelem = str_nelem(val, EGP_MAX, names);
2814 gmx_fatal(FARGS, "The number of groups for %s is odd", option);
2816 nr = groups->grps[egcENER].nr;
2818 for (i = 0; i < nelem/2; i++)
2822 gmx_strcasecmp(names[2*i], *(gnames[groups->grps[egcENER].nm_ind[j]])))
2828 gmx_fatal(FARGS, "%s in %s is not an energy group\n",
2829 names[2*i], option);
2833 gmx_strcasecmp(names[2*i+1], *(gnames[groups->grps[egcENER].nm_ind[k]])))
2839 gmx_fatal(FARGS, "%s in %s is not an energy group\n",
2840 names[2*i+1], option);
2842 if ((j < nr) && (k < nr))
2844 ir->opts.egp_flags[nr*j+k] |= flag;
2845 ir->opts.egp_flags[nr*k+j] |= flag;
2853 void do_index(const char* mdparin, const char *ndx,
2856 t_inputrec *ir, rvec *v,
2860 gmx_groups_t *groups;
2864 char warnbuf[STRLEN], **gnames;
2865 int nr, ntcg, ntau_t, nref_t, nacc, nofg, nSA, nSA_points, nSA_time, nSA_temp;
2868 int nacg, nfreeze, nfrdim, nenergy, nvcm, nuser;
2869 char *ptr1[MAXPTR], *ptr2[MAXPTR], *ptr3[MAXPTR];
2870 int i, j, k, restnm;
2872 gmx_bool bExcl, bTable, bSetTCpar, bAnneal, bRest;
2873 int nQMmethod, nQMbasis, nQMcharge, nQMmult, nbSH, nCASorb, nCASelec,
2874 nSAon, nSAoff, nSAsteps, nQMg, nbOPT, nbTS;
2875 char warn_buf[STRLEN];
2879 fprintf(stderr, "processing index file...\n");
2885 snew(grps->index, 1);
2887 atoms_all = gmx_mtop_global_atoms(mtop);
2888 analyse(&atoms_all, grps, &gnames, FALSE, TRUE);
2889 free_t_atoms(&atoms_all, FALSE);
2893 grps = init_index(ndx, &gnames);
2896 groups = &mtop->groups;
2897 natoms = mtop->natoms;
2898 symtab = &mtop->symtab;
2900 snew(groups->grpname, grps->nr+1);
2902 for (i = 0; (i < grps->nr); i++)
2904 groups->grpname[i] = put_symtab(symtab, gnames[i]);
2906 groups->grpname[i] = put_symtab(symtab, "rest");
2908 srenew(gnames, grps->nr+1);
2909 gnames[restnm] = *(groups->grpname[i]);
2910 groups->ngrpname = grps->nr+1;
2912 set_warning_line(wi, mdparin, -1);
2914 ntau_t = str_nelem(tau_t, MAXPTR, ptr1);
2915 nref_t = str_nelem(ref_t, MAXPTR, ptr2);
2916 ntcg = str_nelem(tcgrps, MAXPTR, ptr3);
2917 if ((ntau_t != ntcg) || (nref_t != ntcg))
2919 gmx_fatal(FARGS, "Invalid T coupling input: %d groups, %d ref-t values and "
2920 "%d tau-t values", ntcg, nref_t, ntau_t);
2923 bSetTCpar = (ir->etc || EI_SD(ir->eI) || ir->eI == eiBD || EI_TPI(ir->eI));
2924 do_numbering(natoms, groups, ntcg, ptr3, grps, gnames, egcTC,
2925 restnm, bSetTCpar ? egrptpALL : egrptpALL_GENREST, bVerbose, wi);
2926 nr = groups->grps[egcTC].nr;
2928 snew(ir->opts.nrdf, nr);
2929 snew(ir->opts.tau_t, nr);
2930 snew(ir->opts.ref_t, nr);
2931 if (ir->eI == eiBD && ir->bd_fric == 0)
2933 fprintf(stderr, "bd-fric=0, so tau-t will be used as the inverse friction constant(s)\n");
2940 gmx_fatal(FARGS, "Not enough ref-t and tau-t values!");
2944 for (i = 0; (i < nr); i++)
2946 ir->opts.tau_t[i] = strtod(ptr1[i], NULL);
2947 if ((ir->eI == eiBD || ir->eI == eiSD2) && ir->opts.tau_t[i] <= 0)
2949 sprintf(warn_buf, "With integrator %s tau-t should be larger than 0", ei_names[ir->eI]);
2950 warning_error(wi, warn_buf);
2953 if (ir->etc != etcVRESCALE && ir->opts.tau_t[i] == 0)
2955 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.");
2958 if (ir->opts.tau_t[i] >= 0)
2960 tau_min = min(tau_min, ir->opts.tau_t[i]);
2963 if (ir->etc != etcNO && ir->nsttcouple == -1)
2965 ir->nsttcouple = ir_optimal_nsttcouple(ir);
2970 if ((ir->etc == etcNOSEHOOVER) && (ir->epc == epcBERENDSEN))
2972 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");
2974 if ((ir->epc == epcMTTK) && (ir->etc > etcNO))
2976 if (ir->nstpcouple != ir->nsttcouple)
2978 int mincouple = min(ir->nstpcouple, ir->nsttcouple);
2979 ir->nstpcouple = ir->nsttcouple = mincouple;
2980 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);
2981 warning_note(wi, warn_buf);
2985 /* velocity verlet with averaged kinetic energy KE = 0.5*(v(t+1/2) - v(t-1/2)) is implemented
2986 primarily for testing purposes, and does not work with temperature coupling other than 1 */
2988 if (ETC_ANDERSEN(ir->etc))
2990 if (ir->nsttcouple != 1)
2993 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");
2994 warning_note(wi, warn_buf);
2997 nstcmin = tcouple_min_integration_steps(ir->etc);
3000 if (tau_min/(ir->delta_t*ir->nsttcouple) < nstcmin)
3002 sprintf(warn_buf, "For proper integration of the %s thermostat, tau-t (%g) should be at least %d times larger than nsttcouple*dt (%g)",
3003 ETCOUPLTYPE(ir->etc),
3005 ir->nsttcouple*ir->delta_t);
3006 warning(wi, warn_buf);
3009 for (i = 0; (i < nr); i++)
3011 ir->opts.ref_t[i] = strtod(ptr2[i], NULL);
3012 if (ir->opts.ref_t[i] < 0)
3014 gmx_fatal(FARGS, "ref-t for group %d negative", i);
3017 /* set the lambda mc temperature to the md integrator temperature (which should be defined
3018 if we are in this conditional) if mc_temp is negative */
3019 if (ir->expandedvals->mc_temp < 0)
3021 ir->expandedvals->mc_temp = ir->opts.ref_t[0]; /*for now, set to the first reft */
3025 /* Simulated annealing for each group. There are nr groups */
3026 nSA = str_nelem(anneal, MAXPTR, ptr1);
3027 if (nSA == 1 && (ptr1[0][0] == 'n' || ptr1[0][0] == 'N'))
3031 if (nSA > 0 && nSA != nr)
3033 gmx_fatal(FARGS, "Not enough annealing values: %d (for %d groups)\n", nSA, nr);
3037 snew(ir->opts.annealing, nr);
3038 snew(ir->opts.anneal_npoints, nr);
3039 snew(ir->opts.anneal_time, nr);
3040 snew(ir->opts.anneal_temp, nr);
3041 for (i = 0; i < nr; i++)
3043 ir->opts.annealing[i] = eannNO;
3044 ir->opts.anneal_npoints[i] = 0;
3045 ir->opts.anneal_time[i] = NULL;
3046 ir->opts.anneal_temp[i] = NULL;
3051 for (i = 0; i < nr; i++)
3053 if (ptr1[i][0] == 'n' || ptr1[i][0] == 'N')
3055 ir->opts.annealing[i] = eannNO;
3057 else if (ptr1[i][0] == 's' || ptr1[i][0] == 'S')
3059 ir->opts.annealing[i] = eannSINGLE;
3062 else if (ptr1[i][0] == 'p' || ptr1[i][0] == 'P')
3064 ir->opts.annealing[i] = eannPERIODIC;
3070 /* Read the other fields too */
3071 nSA_points = str_nelem(anneal_npoints, MAXPTR, ptr1);
3072 if (nSA_points != nSA)
3074 gmx_fatal(FARGS, "Found %d annealing-npoints values for %d groups\n", nSA_points, nSA);
3076 for (k = 0, i = 0; i < nr; i++)
3078 ir->opts.anneal_npoints[i] = strtol(ptr1[i], NULL, 10);
3079 if (ir->opts.anneal_npoints[i] == 1)
3081 gmx_fatal(FARGS, "Please specify at least a start and an end point for annealing\n");
3083 snew(ir->opts.anneal_time[i], ir->opts.anneal_npoints[i]);
3084 snew(ir->opts.anneal_temp[i], ir->opts.anneal_npoints[i]);
3085 k += ir->opts.anneal_npoints[i];
3088 nSA_time = str_nelem(anneal_time, MAXPTR, ptr1);
3091 gmx_fatal(FARGS, "Found %d annealing-time values, wanter %d\n", nSA_time, k);
3093 nSA_temp = str_nelem(anneal_temp, MAXPTR, ptr2);
3096 gmx_fatal(FARGS, "Found %d annealing-temp values, wanted %d\n", nSA_temp, k);
3099 for (i = 0, k = 0; i < nr; i++)
3102 for (j = 0; j < ir->opts.anneal_npoints[i]; j++)
3104 ir->opts.anneal_time[i][j] = strtod(ptr1[k], NULL);
3105 ir->opts.anneal_temp[i][j] = strtod(ptr2[k], NULL);
3108 if (ir->opts.anneal_time[i][0] > (ir->init_t+GMX_REAL_EPS))
3110 gmx_fatal(FARGS, "First time point for annealing > init_t.\n");
3116 if (ir->opts.anneal_time[i][j] < ir->opts.anneal_time[i][j-1])
3118 gmx_fatal(FARGS, "Annealing timepoints out of order: t=%f comes after t=%f\n",
3119 ir->opts.anneal_time[i][j], ir->opts.anneal_time[i][j-1]);
3122 if (ir->opts.anneal_temp[i][j] < 0)
3124 gmx_fatal(FARGS, "Found negative temperature in annealing: %f\n", ir->opts.anneal_temp[i][j]);
3129 /* Print out some summary information, to make sure we got it right */
3130 for (i = 0, k = 0; i < nr; i++)
3132 if (ir->opts.annealing[i] != eannNO)
3134 j = groups->grps[egcTC].nm_ind[i];
3135 fprintf(stderr, "Simulated annealing for group %s: %s, %d timepoints\n",
3136 *(groups->grpname[j]), eann_names[ir->opts.annealing[i]],
3137 ir->opts.anneal_npoints[i]);
3138 fprintf(stderr, "Time (ps) Temperature (K)\n");
3139 /* All terms except the last one */
3140 for (j = 0; j < (ir->opts.anneal_npoints[i]-1); j++)
3142 fprintf(stderr, "%9.1f %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3145 /* Finally the last one */
3146 j = ir->opts.anneal_npoints[i]-1;
3147 if (ir->opts.annealing[i] == eannSINGLE)
3149 fprintf(stderr, "%9.1f- %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3153 fprintf(stderr, "%9.1f %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3154 if (fabs(ir->opts.anneal_temp[i][j]-ir->opts.anneal_temp[i][0]) > GMX_REAL_EPS)
3156 warning_note(wi, "There is a temperature jump when your annealing loops back.\n");
3165 if (ir->ePull != epullNO)
3167 make_pull_groups(ir->pull, pull_grp, grps, gnames);
3169 make_pull_coords(ir->pull);
3174 make_rotation_groups(ir->rot, rot_grp, grps, gnames);
3177 nacc = str_nelem(acc, MAXPTR, ptr1);
3178 nacg = str_nelem(accgrps, MAXPTR, ptr2);
3179 if (nacg*DIM != nacc)
3181 gmx_fatal(FARGS, "Invalid Acceleration input: %d groups and %d acc. values",
3184 do_numbering(natoms, groups, nacg, ptr2, grps, gnames, egcACC,
3185 restnm, egrptpALL_GENREST, bVerbose, wi);
3186 nr = groups->grps[egcACC].nr;
3187 snew(ir->opts.acc, nr);
3188 ir->opts.ngacc = nr;
3190 for (i = k = 0; (i < nacg); i++)
3192 for (j = 0; (j < DIM); j++, k++)
3194 ir->opts.acc[i][j] = strtod(ptr1[k], NULL);
3197 for (; (i < nr); i++)
3199 for (j = 0; (j < DIM); j++)
3201 ir->opts.acc[i][j] = 0;
3205 nfrdim = str_nelem(frdim, MAXPTR, ptr1);
3206 nfreeze = str_nelem(freeze, MAXPTR, ptr2);
3207 if (nfrdim != DIM*nfreeze)
3209 gmx_fatal(FARGS, "Invalid Freezing input: %d groups and %d freeze values",
3212 do_numbering(natoms, groups, nfreeze, ptr2, grps, gnames, egcFREEZE,
3213 restnm, egrptpALL_GENREST, bVerbose, wi);
3214 nr = groups->grps[egcFREEZE].nr;
3215 ir->opts.ngfrz = nr;
3216 snew(ir->opts.nFreeze, nr);
3217 for (i = k = 0; (i < nfreeze); i++)
3219 for (j = 0; (j < DIM); j++, k++)
3221 ir->opts.nFreeze[i][j] = (gmx_strncasecmp(ptr1[k], "Y", 1) == 0);
3222 if (!ir->opts.nFreeze[i][j])
3224 if (gmx_strncasecmp(ptr1[k], "N", 1) != 0)
3226 sprintf(warnbuf, "Please use Y(ES) or N(O) for freezedim only "
3227 "(not %s)", ptr1[k]);
3228 warning(wi, warn_buf);
3233 for (; (i < nr); i++)
3235 for (j = 0; (j < DIM); j++)
3237 ir->opts.nFreeze[i][j] = 0;
3241 nenergy = str_nelem(energy, MAXPTR, ptr1);
3242 do_numbering(natoms, groups, nenergy, ptr1, grps, gnames, egcENER,
3243 restnm, egrptpALL_GENREST, bVerbose, wi);
3244 add_wall_energrps(groups, ir->nwall, symtab);
3245 ir->opts.ngener = groups->grps[egcENER].nr;
3246 nvcm = str_nelem(vcm, MAXPTR, ptr1);
3248 do_numbering(natoms, groups, nvcm, ptr1, grps, gnames, egcVCM,
3249 restnm, nvcm == 0 ? egrptpALL_GENREST : egrptpPART, bVerbose, wi);
3252 warning(wi, "Some atoms are not part of any center of mass motion removal group.\n"
3253 "This may lead to artifacts.\n"
3254 "In most cases one should use one group for the whole system.");
3257 /* Now we have filled the freeze struct, so we can calculate NRDF */
3258 calc_nrdf(mtop, ir, gnames);
3264 /* Must check per group! */
3265 for (i = 0; (i < ir->opts.ngtc); i++)
3267 ntot += ir->opts.nrdf[i];
3269 if (ntot != (DIM*natoms))
3271 fac = sqrt(ntot/(DIM*natoms));
3274 fprintf(stderr, "Scaling velocities by a factor of %.3f to account for constraints\n"
3275 "and removal of center of mass motion\n", fac);
3277 for (i = 0; (i < natoms); i++)
3279 svmul(fac, v[i], v[i]);
3284 nuser = str_nelem(user1, MAXPTR, ptr1);
3285 do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcUser1,
3286 restnm, egrptpALL_GENREST, bVerbose, wi);
3287 nuser = str_nelem(user2, MAXPTR, ptr1);
3288 do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcUser2,
3289 restnm, egrptpALL_GENREST, bVerbose, wi);
3290 nuser = str_nelem(xtc_grps, MAXPTR, ptr1);
3291 do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcXTC,
3292 restnm, egrptpONE, bVerbose, wi);
3293 nofg = str_nelem(orirefitgrp, MAXPTR, ptr1);
3294 do_numbering(natoms, groups, nofg, ptr1, grps, gnames, egcORFIT,
3295 restnm, egrptpALL_GENREST, bVerbose, wi);
3297 /* QMMM input processing */
3298 nQMg = str_nelem(QMMM, MAXPTR, ptr1);
3299 nQMmethod = str_nelem(QMmethod, MAXPTR, ptr2);
3300 nQMbasis = str_nelem(QMbasis, MAXPTR, ptr3);
3301 if ((nQMmethod != nQMg) || (nQMbasis != nQMg))
3303 gmx_fatal(FARGS, "Invalid QMMM input: %d groups %d basissets"
3304 " and %d methods\n", nQMg, nQMbasis, nQMmethod);
3306 /* group rest, if any, is always MM! */
3307 do_numbering(natoms, groups, nQMg, ptr1, grps, gnames, egcQMMM,
3308 restnm, egrptpALL_GENREST, bVerbose, wi);
3309 nr = nQMg; /*atoms->grps[egcQMMM].nr;*/
3310 ir->opts.ngQM = nQMg;
3311 snew(ir->opts.QMmethod, nr);
3312 snew(ir->opts.QMbasis, nr);
3313 for (i = 0; i < nr; i++)
3315 /* input consists of strings: RHF CASSCF PM3 .. These need to be
3316 * converted to the corresponding enum in names.c
3318 ir->opts.QMmethod[i] = search_QMstring(ptr2[i], eQMmethodNR,
3320 ir->opts.QMbasis[i] = search_QMstring(ptr3[i], eQMbasisNR,
3324 nQMmult = str_nelem(QMmult, MAXPTR, ptr1);
3325 nQMcharge = str_nelem(QMcharge, MAXPTR, ptr2);
3326 nbSH = str_nelem(bSH, MAXPTR, ptr3);
3327 snew(ir->opts.QMmult, nr);
3328 snew(ir->opts.QMcharge, nr);
3329 snew(ir->opts.bSH, nr);
3331 for (i = 0; i < nr; i++)
3333 ir->opts.QMmult[i] = strtol(ptr1[i], NULL, 10);
3334 ir->opts.QMcharge[i] = strtol(ptr2[i], NULL, 10);
3335 ir->opts.bSH[i] = (gmx_strncasecmp(ptr3[i], "Y", 1) == 0);
3338 nCASelec = str_nelem(CASelectrons, MAXPTR, ptr1);
3339 nCASorb = str_nelem(CASorbitals, MAXPTR, ptr2);
3340 snew(ir->opts.CASelectrons, nr);
3341 snew(ir->opts.CASorbitals, nr);
3342 for (i = 0; i < nr; i++)
3344 ir->opts.CASelectrons[i] = strtol(ptr1[i], NULL, 10);
3345 ir->opts.CASorbitals[i] = strtol(ptr2[i], NULL, 10);
3347 /* special optimization options */
3349 nbOPT = str_nelem(bOPT, MAXPTR, ptr1);
3350 nbTS = str_nelem(bTS, MAXPTR, ptr2);
3351 snew(ir->opts.bOPT, nr);
3352 snew(ir->opts.bTS, nr);
3353 for (i = 0; i < nr; i++)
3355 ir->opts.bOPT[i] = (gmx_strncasecmp(ptr1[i], "Y", 1) == 0);
3356 ir->opts.bTS[i] = (gmx_strncasecmp(ptr2[i], "Y", 1) == 0);
3358 nSAon = str_nelem(SAon, MAXPTR, ptr1);
3359 nSAoff = str_nelem(SAoff, MAXPTR, ptr2);
3360 nSAsteps = str_nelem(SAsteps, MAXPTR, ptr3);
3361 snew(ir->opts.SAon, nr);
3362 snew(ir->opts.SAoff, nr);
3363 snew(ir->opts.SAsteps, nr);
3365 for (i = 0; i < nr; i++)
3367 ir->opts.SAon[i] = strtod(ptr1[i], NULL);
3368 ir->opts.SAoff[i] = strtod(ptr2[i], NULL);
3369 ir->opts.SAsteps[i] = strtol(ptr3[i], NULL, 10);
3371 /* end of QMMM input */
3375 for (i = 0; (i < egcNR); i++)
3377 fprintf(stderr, "%-16s has %d element(s):", gtypes[i], groups->grps[i].nr);
3378 for (j = 0; (j < groups->grps[i].nr); j++)
3380 fprintf(stderr, " %s", *(groups->grpname[groups->grps[i].nm_ind[j]]));
3382 fprintf(stderr, "\n");
3386 nr = groups->grps[egcENER].nr;
3387 snew(ir->opts.egp_flags, nr*nr);
3389 bExcl = do_egp_flag(ir, groups, "energygrp-excl", egpexcl, EGP_EXCL);
3390 if (bExcl && ir->cutoff_scheme == ecutsVERLET)
3392 warning_error(wi, "Energy group exclusions are not (yet) implemented for the Verlet scheme");
3394 if (bExcl && EEL_FULL(ir->coulombtype))
3396 warning(wi, "Can not exclude the lattice Coulomb energy between energy groups");
3399 bTable = do_egp_flag(ir, groups, "energygrp-table", egptable, EGP_TABLE);
3400 if (bTable && !(ir->vdwtype == evdwUSER) &&
3401 !(ir->coulombtype == eelUSER) && !(ir->coulombtype == eelPMEUSER) &&
3402 !(ir->coulombtype == eelPMEUSERSWITCH))
3404 gmx_fatal(FARGS, "Can only have energy group pair tables in combination with user tables for VdW and/or Coulomb");
3407 decode_cos(efield_x, &(ir->ex[XX]));
3408 decode_cos(efield_xt, &(ir->et[XX]));
3409 decode_cos(efield_y, &(ir->ex[YY]));
3410 decode_cos(efield_yt, &(ir->et[YY]));
3411 decode_cos(efield_z, &(ir->ex[ZZ]));
3412 decode_cos(efield_zt, &(ir->et[ZZ]));
3416 do_adress_index(ir->adress, groups, gnames, &(ir->opts), wi);
3419 for (i = 0; (i < grps->nr); i++)
3431 static void check_disre(gmx_mtop_t *mtop)
3433 gmx_ffparams_t *ffparams;
3434 t_functype *functype;
3436 int i, ndouble, ftype;
3437 int label, old_label;
3439 if (gmx_mtop_ftype_count(mtop, F_DISRES) > 0)
3441 ffparams = &mtop->ffparams;
3442 functype = ffparams->functype;
3443 ip = ffparams->iparams;
3446 for (i = 0; i < ffparams->ntypes; i++)
3448 ftype = functype[i];
3449 if (ftype == F_DISRES)
3451 label = ip[i].disres.label;
3452 if (label == old_label)
3454 fprintf(stderr, "Distance restraint index %d occurs twice\n", label);
3462 gmx_fatal(FARGS, "Found %d double distance restraint indices,\n"
3463 "probably the parameters for multiple pairs in one restraint "
3464 "are not identical\n", ndouble);
3469 static gmx_bool absolute_reference(t_inputrec *ir, gmx_mtop_t *sys,
3470 gmx_bool posres_only,
3474 gmx_mtop_ilistloop_t iloop;
3484 for (d = 0; d < DIM; d++)
3486 AbsRef[d] = (d < ndof_com(ir) ? 0 : 1);
3488 /* Check for freeze groups */
3489 for (g = 0; g < ir->opts.ngfrz; g++)
3491 for (d = 0; d < DIM; d++)
3493 if (ir->opts.nFreeze[g][d] != 0)
3501 /* Check for position restraints */
3502 iloop = gmx_mtop_ilistloop_init(sys);
3503 while (gmx_mtop_ilistloop_next(iloop, &ilist, &nmol))
3506 (AbsRef[XX] == 0 || AbsRef[YY] == 0 || AbsRef[ZZ] == 0))
3508 for (i = 0; i < ilist[F_POSRES].nr; i += 2)
3510 pr = &sys->ffparams.iparams[ilist[F_POSRES].iatoms[i]];
3511 for (d = 0; d < DIM; d++)
3513 if (pr->posres.fcA[d] != 0)
3519 for (i = 0; i < ilist[F_FBPOSRES].nr; i += 2)
3521 /* Check for flat-bottom posres */
3522 pr = &sys->ffparams.iparams[ilist[F_FBPOSRES].iatoms[i]];
3523 if (pr->fbposres.k != 0)
3525 switch (pr->fbposres.geom)
3527 case efbposresSPHERE:
3528 AbsRef[XX] = AbsRef[YY] = AbsRef[ZZ] = 1;
3530 case efbposresCYLINDER:
3531 AbsRef[XX] = AbsRef[YY] = 1;
3533 case efbposresX: /* d=XX */
3534 case efbposresY: /* d=YY */
3535 case efbposresZ: /* d=ZZ */
3536 d = pr->fbposres.geom - efbposresX;
3540 gmx_fatal(FARGS, " Invalid geometry for flat-bottom position restraint.\n"
3541 "Expected nr between 1 and %d. Found %d\n", efbposresNR-1,
3549 return (AbsRef[XX] != 0 && AbsRef[YY] != 0 && AbsRef[ZZ] != 0);
3553 check_combination_rule_differences(const gmx_mtop_t *mtop, int state,
3554 gmx_bool *bC6ParametersWorkWithGeometricRules,
3555 gmx_bool *bC6ParametersWorkWithLBRules,
3556 gmx_bool *bLBRulesPossible)
3558 int ntypes, tpi, tpj, thisLBdiff, thisgeomdiff;
3561 #if (defined SIZEOF_LONG_LONG_INT) && (SIZEOF_LONG_LONG_INT >= 8)
3562 long long int npair, npair_ij;
3564 double npair, npair_ij;
3566 double geometricdiff, LBdiff;
3567 double c6i, c6j, c12i, c12j;
3568 double c6, c6_geometric, c6_LB;
3569 double sigmai, sigmaj, epsi, epsj;
3570 gmx_bool bCanDoLBRules, bCanDoGeometricRules;
3573 /* A tolerance of 1e-5 seems reasonable for (possibly hand-typed)
3574 * force-field floating point parameters.
3577 ptr = getenv("GMX_LJCOMB_TOL");
3582 sscanf(ptr, "%lf", &dbl);
3586 *bC6ParametersWorkWithLBRules = TRUE;
3587 *bC6ParametersWorkWithGeometricRules = TRUE;
3588 bCanDoLBRules = TRUE;
3589 bCanDoGeometricRules = TRUE;
3591 ntypes = mtop->ffparams.atnr;
3592 snew(typecount, ntypes);
3593 gmx_mtop_count_atomtypes(mtop, state, typecount);
3594 geometricdiff = LBdiff = 0.0;
3595 *bLBRulesPossible = TRUE;
3596 for (tpi = 0; tpi < ntypes; ++tpi)
3598 c6i = mtop->ffparams.iparams[(ntypes + 1) * tpi].lj.c6;
3599 c12i = mtop->ffparams.iparams[(ntypes + 1) * tpi].lj.c12;
3600 for (tpj = tpi; tpj < ntypes; ++tpj)
3602 c6j = mtop->ffparams.iparams[(ntypes + 1) * tpj].lj.c6;
3603 c12j = mtop->ffparams.iparams[(ntypes + 1) * tpj].lj.c12;
3604 c6 = mtop->ffparams.iparams[ntypes * tpi + tpj].lj.c6;
3605 c6_geometric = sqrt(c6i * c6j);
3606 if (!gmx_numzero(c6_geometric))
3608 if (!gmx_numzero(c12i) && !gmx_numzero(c12j))
3610 sigmai = pow(c12i / c6i, 1.0/6.0);
3611 sigmaj = pow(c12j / c6j, 1.0/6.0);
3612 epsi = c6i * c6i /(4.0 * c12i);
3613 epsj = c6j * c6j /(4.0 * c12j);
3614 c6_LB = 4.0 * pow(epsi * epsj, 1.0/2.0) * pow(0.5 * (sigmai + sigmaj), 6);
3618 *bLBRulesPossible = FALSE;
3619 c6_LB = c6_geometric;
3621 bCanDoLBRules = gmx_within_tol(c6_LB, c6, tol);
3624 if (FALSE == bCanDoLBRules)
3626 *bC6ParametersWorkWithLBRules = FALSE;
3629 bCanDoGeometricRules = gmx_within_tol(c6_geometric, c6, tol);
3631 if (FALSE == bCanDoGeometricRules)
3633 *bC6ParametersWorkWithGeometricRules = FALSE;
3641 check_combination_rules(const t_inputrec *ir, const gmx_mtop_t *mtop,
3645 gmx_bool bLBRulesPossible, bC6ParametersWorkWithGeometricRules, bC6ParametersWorkWithLBRules;
3647 check_combination_rule_differences(mtop, 0,
3648 &bC6ParametersWorkWithGeometricRules,
3649 &bC6ParametersWorkWithLBRules,
3651 if (ir->ljpme_combination_rule == eljpmeLB)
3653 if (FALSE == bC6ParametersWorkWithLBRules || FALSE == bLBRulesPossible)
3655 warning(wi, "You are using arithmetic-geometric combination rules "
3656 "in LJ-PME, but your non-bonded C6 parameters do not "
3657 "follow these rules.");
3662 if (FALSE == bC6ParametersWorkWithGeometricRules)
3664 if (ir->eDispCorr != edispcNO)
3666 warning_note(wi, "You are using geometric combination rules in "
3667 "LJ-PME, but your non-bonded C6 parameters do "
3668 "not follow these rules. "
3669 "If your force field uses Lorentz-Berthelot combination rules this "
3670 "will introduce small errors in the forces and energies in "
3671 "your simulations. Dispersion correction will correct total "
3672 "energy and/or pressure, but not forces or surface tensions. "
3673 "Please check the LJ-PME section in the manual "
3674 "before proceeding further.");
3678 warning_note(wi, "You are using geometric combination rules in "
3679 "LJ-PME, but your non-bonded C6 parameters do "
3680 "not follow these rules. "
3681 "If your force field uses Lorentz-Berthelot combination rules this "
3682 "will introduce small errors in the forces and energies in "
3683 "your simulations. Consider using dispersion correction "
3684 "for the total energy and pressure. "
3685 "Please check the LJ-PME section in the manual "
3686 "before proceeding further.");
3692 void triple_check(const char *mdparin, t_inputrec *ir, gmx_mtop_t *sys,
3696 int i, m, c, nmol, npct;
3697 gmx_bool bCharge, bAcc;
3698 real gdt_max, *mgrp, mt;
3700 gmx_mtop_atomloop_block_t aloopb;
3701 gmx_mtop_atomloop_all_t aloop;
3704 char warn_buf[STRLEN];
3706 set_warning_line(wi, mdparin, -1);
3708 if (EI_DYNAMICS(ir->eI) && !EI_SD(ir->eI) && ir->eI != eiBD &&
3709 ir->comm_mode == ecmNO &&
3710 !(absolute_reference(ir, sys, FALSE, AbsRef) || ir->nsteps <= 10))
3712 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");
3715 /* Check for pressure coupling with absolute position restraints */
3716 if (ir->epc != epcNO && ir->refcoord_scaling == erscNO)
3718 absolute_reference(ir, sys, TRUE, AbsRef);
3720 for (m = 0; m < DIM; m++)
3722 if (AbsRef[m] && norm2(ir->compress[m]) > 0)
3724 warning(wi, "You are using pressure coupling with absolute position restraints, this will give artifacts. Use the refcoord_scaling option.");
3732 aloopb = gmx_mtop_atomloop_block_init(sys);
3733 while (gmx_mtop_atomloop_block_next(aloopb, &atom, &nmol))
3735 if (atom->q != 0 || atom->qB != 0)
3743 if (EEL_FULL(ir->coulombtype))
3746 "You are using full electrostatics treatment %s for a system without charges.\n"
3747 "This costs a lot of performance for just processing zeros, consider using %s instead.\n",
3748 EELTYPE(ir->coulombtype), EELTYPE(eelCUT));
3749 warning(wi, err_buf);
3754 if (ir->coulombtype == eelCUT && ir->rcoulomb > 0 && !ir->implicit_solvent)
3757 "You are using a plain Coulomb cut-off, which might produce artifacts.\n"
3758 "You might want to consider using %s electrostatics.\n",
3760 warning_note(wi, err_buf);
3764 /* Check if combination rules used in LJ-PME are the same as in the force field */
3765 if (EVDW_PME(ir->vdwtype))
3767 check_combination_rules(ir, sys, wi);
3770 /* Generalized reaction field */
3771 if (ir->opts.ngtc == 0)
3773 sprintf(err_buf, "No temperature coupling while using coulombtype %s",
3775 CHECK(ir->coulombtype == eelGRF);
3779 sprintf(err_buf, "When using coulombtype = %s"
3780 " ref-t for temperature coupling should be > 0",
3782 CHECK((ir->coulombtype == eelGRF) && (ir->opts.ref_t[0] <= 0));
3785 if (ir->eI == eiSD1 &&
3786 (gmx_mtop_ftype_count(sys, F_CONSTR) > 0 ||
3787 gmx_mtop_ftype_count(sys, F_SETTLE) > 0))
3789 sprintf(warn_buf, "With constraints integrator %s is less accurate, consider using %s instead", ei_names[ir->eI], ei_names[eiSD2]);
3790 warning_note(wi, warn_buf);
3794 for (i = 0; (i < sys->groups.grps[egcACC].nr); i++)
3796 for (m = 0; (m < DIM); m++)
3798 if (fabs(ir->opts.acc[i][m]) > 1e-6)
3807 snew(mgrp, sys->groups.grps[egcACC].nr);
3808 aloop = gmx_mtop_atomloop_all_init(sys);
3809 while (gmx_mtop_atomloop_all_next(aloop, &i, &atom))
3811 mgrp[ggrpnr(&sys->groups, egcACC, i)] += atom->m;
3814 for (i = 0; (i < sys->groups.grps[egcACC].nr); i++)
3816 for (m = 0; (m < DIM); m++)
3818 acc[m] += ir->opts.acc[i][m]*mgrp[i];
3822 for (m = 0; (m < DIM); m++)
3824 if (fabs(acc[m]) > 1e-6)
3826 const char *dim[DIM] = { "X", "Y", "Z" };
3828 "Net Acceleration in %s direction, will %s be corrected\n",
3829 dim[m], ir->nstcomm != 0 ? "" : "not");
3830 if (ir->nstcomm != 0 && m < ndof_com(ir))
3833 for (i = 0; (i < sys->groups.grps[egcACC].nr); i++)
3835 ir->opts.acc[i][m] -= acc[m];
3843 if (ir->efep != efepNO && ir->fepvals->sc_alpha != 0 &&
3844 !gmx_within_tol(sys->ffparams.reppow, 12.0, 10*GMX_DOUBLE_EPS))
3846 gmx_fatal(FARGS, "Soft-core interactions are only supported with VdW repulsion power 12");
3849 if (ir->ePull != epullNO)
3851 gmx_bool bPullAbsoluteRef;
3853 bPullAbsoluteRef = FALSE;
3854 for (i = 0; i < ir->pull->ncoord; i++)
3856 bPullAbsoluteRef = bPullAbsoluteRef ||
3857 ir->pull->coord[i].group[0] == 0 ||
3858 ir->pull->coord[i].group[1] == 0;
3860 if (bPullAbsoluteRef)
3862 absolute_reference(ir, sys, FALSE, AbsRef);
3863 for (m = 0; m < DIM; m++)
3865 if (ir->pull->dim[m] && !AbsRef[m])
3867 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.");
3873 if (ir->pull->eGeom == epullgDIRPBC)
3875 for (i = 0; i < 3; i++)
3877 for (m = 0; m <= i; m++)
3879 if ((ir->epc != epcNO && ir->compress[i][m] != 0) ||
3880 ir->deform[i][m] != 0)
3882 for (c = 0; c < ir->pull->ncoord; c++)
3884 if (ir->pull->coord[c].vec[m] != 0)
3886 gmx_fatal(FARGS, "Can not have dynamic box while using pull geometry '%s' (dim %c)", EPULLGEOM(ir->pull->eGeom), 'x'+m);
3898 void double_check(t_inputrec *ir, matrix box, gmx_bool bConstr, warninp_t wi)
3902 char warn_buf[STRLEN];
3905 ptr = check_box(ir->ePBC, box);
3908 warning_error(wi, ptr);
3911 if (bConstr && ir->eConstrAlg == econtSHAKE)
3913 if (ir->shake_tol <= 0.0)
3915 sprintf(warn_buf, "ERROR: shake-tol must be > 0 instead of %g\n",
3917 warning_error(wi, warn_buf);
3920 if (IR_TWINRANGE(*ir) && ir->nstlist > 1)
3922 sprintf(warn_buf, "With twin-range cut-off's and SHAKE the virial and the pressure are incorrect.");
3923 if (ir->epc == epcNO)
3925 warning(wi, warn_buf);
3929 warning_error(wi, warn_buf);
3934 if ( (ir->eConstrAlg == econtLINCS) && bConstr)
3936 /* If we have Lincs constraints: */
3937 if (ir->eI == eiMD && ir->etc == etcNO &&
3938 ir->eConstrAlg == econtLINCS && ir->nLincsIter == 1)
3940 sprintf(warn_buf, "For energy conservation with LINCS, lincs_iter should be 2 or larger.\n");
3941 warning_note(wi, warn_buf);
3944 if ((ir->eI == eiCG || ir->eI == eiLBFGS) && (ir->nProjOrder < 8))
3946 sprintf(warn_buf, "For accurate %s with LINCS constraints, lincs-order should be 8 or more.", ei_names[ir->eI]);
3947 warning_note(wi, warn_buf);
3949 if (ir->epc == epcMTTK)
3951 warning_error(wi, "MTTK not compatible with lincs -- use shake instead.");
3955 if (ir->LincsWarnAngle > 90.0)
3957 sprintf(warn_buf, "lincs-warnangle can not be larger than 90 degrees, setting it to 90.\n");
3958 warning(wi, warn_buf);
3959 ir->LincsWarnAngle = 90.0;
3962 if (ir->ePBC != epbcNONE)
3964 if (ir->nstlist == 0)
3966 warning(wi, "With nstlist=0 atoms are only put into the box at step 0, therefore drifting atoms might cause the simulation to crash.");
3968 bTWIN = (ir->rlistlong > ir->rlist);
3969 if (ir->ns_type == ensGRID)
3971 if (sqr(ir->rlistlong) >= max_cutoff2(ir->ePBC, box))
3973 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",
3974 bTWIN ? (ir->rcoulomb == ir->rlistlong ? "rcoulomb" : "rvdw") : "rlist");
3975 warning_error(wi, warn_buf);
3980 min_size = min(box[XX][XX], min(box[YY][YY], box[ZZ][ZZ]));
3981 if (2*ir->rlistlong >= min_size)
3983 sprintf(warn_buf, "ERROR: One of the box lengths is smaller than twice the cut-off length. Increase the box size or decrease rlist.");
3984 warning_error(wi, warn_buf);
3987 fprintf(stderr, "Grid search might allow larger cut-off's than simple search with triclinic boxes.");
3994 void check_chargegroup_radii(const gmx_mtop_t *mtop, const t_inputrec *ir,
3998 real rvdw1, rvdw2, rcoul1, rcoul2;
3999 char warn_buf[STRLEN];
4001 calc_chargegroup_radii(mtop, x, &rvdw1, &rvdw2, &rcoul1, &rcoul2);
4005 printf("Largest charge group radii for Van der Waals: %5.3f, %5.3f nm\n",
4010 printf("Largest charge group radii for Coulomb: %5.3f, %5.3f nm\n",
4016 if (rvdw1 + rvdw2 > ir->rlist ||
4017 rcoul1 + rcoul2 > ir->rlist)
4020 "The sum of the two largest charge group radii (%f) "
4021 "is larger than rlist (%f)\n",
4022 max(rvdw1+rvdw2, rcoul1+rcoul2), ir->rlist);
4023 warning(wi, warn_buf);
4027 /* Here we do not use the zero at cut-off macro,
4028 * since user defined interactions might purposely
4029 * not be zero at the cut-off.
4031 if ((EVDW_IS_ZERO_AT_CUTOFF(ir->vdwtype) ||
4032 ir->vdw_modifier != eintmodNONE) &&
4033 rvdw1 + rvdw2 > ir->rlistlong - ir->rvdw)
4035 sprintf(warn_buf, "The sum of the two largest charge group "
4036 "radii (%f) is larger than %s (%f) - rvdw (%f).\n"
4037 "With exact cut-offs, better performance can be "
4038 "obtained with cutoff-scheme = %s, because it "
4039 "does not use charge groups at all.",
4041 ir->rlistlong > ir->rlist ? "rlistlong" : "rlist",
4042 ir->rlistlong, ir->rvdw,
4043 ecutscheme_names[ecutsVERLET]);
4046 warning(wi, warn_buf);
4050 warning_note(wi, warn_buf);
4053 if ((EEL_IS_ZERO_AT_CUTOFF(ir->coulombtype) ||
4054 ir->coulomb_modifier != eintmodNONE) &&
4055 rcoul1 + rcoul2 > ir->rlistlong - ir->rcoulomb)
4057 sprintf(warn_buf, "The sum of the two largest charge group radii (%f) is larger than %s (%f) - rcoulomb (%f).\n"
4058 "With exact cut-offs, better performance can be obtained with cutoff-scheme = %s, because it does not use charge groups at all.",
4060 ir->rlistlong > ir->rlist ? "rlistlong" : "rlist",
4061 ir->rlistlong, ir->rcoulomb,
4062 ecutscheme_names[ecutsVERLET]);
4065 warning(wi, warn_buf);
4069 warning_note(wi, warn_buf);