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 (group: using charge groups, Verlet: particle based cut-offs)");
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(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 */
2330 int search_string(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, nc;
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.
2644 if (pull->eGeom == epullgPOS)
2647 for (i = 0; i < DIM; i++)
2659 for (i = 0; i < pull->ngrp; i++)
2662 if (pull->grp[0].nat > 0)
2664 /* Subtract 1/2 dof from the reference group */
2665 ai = pull->grp[0].ind[0];
2666 if (nrdf_tc[ggrpnr(groups, egcTC, ai)] > 1)
2668 nrdf_tc [ggrpnr(groups, egcTC, ai)] -= 0.5;
2669 nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5;
2673 /* Subtract 1/2 dof from the pulled group */
2674 ai = pull->grp[1+i].ind[0];
2675 nrdf_tc [ggrpnr(groups, egcTC, ai)] -= 0.5*imin;
2676 nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5*imin;
2677 if (nrdf_tc[ggrpnr(groups, egcTC, ai)] < 0)
2679 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)]]);
2684 if (ir->nstcomm != 0)
2686 /* Subtract 3 from the number of degrees of freedom in each vcm group
2687 * when com translation is removed and 6 when rotation is removed
2690 switch (ir->comm_mode)
2693 n_sub = ndof_com(ir);
2700 gmx_incons("Checking comm_mode");
2703 for (i = 0; i < groups->grps[egcTC].nr; i++)
2705 /* Count the number of atoms of TC group i for every VCM group */
2706 for (j = 0; j < groups->grps[egcVCM].nr+1; j++)
2711 for (ai = 0; ai < natoms; ai++)
2713 if (ggrpnr(groups, egcTC, ai) == i)
2715 na_vcm[ggrpnr(groups, egcVCM, ai)]++;
2719 /* Correct for VCM removal according to the fraction of each VCM
2720 * group present in this TC group.
2722 nrdf_uc = nrdf_tc[i];
2725 fprintf(debug, "T-group[%d] nrdf_uc = %g, n_sub = %g\n",
2729 for (j = 0; j < groups->grps[egcVCM].nr+1; j++)
2731 if (nrdf_vcm[j] > n_sub)
2733 nrdf_tc[i] += nrdf_uc*((double)na_vcm[j]/(double)na_tot)*
2734 (nrdf_vcm[j] - n_sub)/nrdf_vcm[j];
2738 fprintf(debug, " nrdf_vcm[%d] = %g, nrdf = %g\n",
2739 j, nrdf_vcm[j], nrdf_tc[i]);
2744 for (i = 0; (i < groups->grps[egcTC].nr); i++)
2746 opts->nrdf[i] = nrdf_tc[i];
2747 if (opts->nrdf[i] < 0)
2752 "Number of degrees of freedom in T-Coupling group %s is %.2f\n",
2753 gnames[groups->grps[egcTC].nm_ind[i]], opts->nrdf[i]);
2762 static void decode_cos(char *s, t_cosines *cosine)
2765 char format[STRLEN], f1[STRLEN];
2777 sscanf(t, "%d", &(cosine->n));
2784 snew(cosine->a, cosine->n);
2785 snew(cosine->phi, cosine->n);
2787 sprintf(format, "%%*d");
2788 for (i = 0; (i < cosine->n); i++)
2791 strcat(f1, "%lf%lf");
2792 if (sscanf(t, f1, &a, &phi) < 2)
2794 gmx_fatal(FARGS, "Invalid input for electric field shift: '%s'", t);
2797 cosine->phi[i] = phi;
2798 strcat(format, "%*lf%*lf");
2805 static gmx_bool do_egp_flag(t_inputrec *ir, gmx_groups_t *groups,
2806 const char *option, const char *val, int flag)
2808 /* The maximum number of energy group pairs would be MAXPTR*(MAXPTR+1)/2.
2809 * But since this is much larger than STRLEN, such a line can not be parsed.
2810 * The real maximum is the number of names that fit in a string: STRLEN/2.
2812 #define EGP_MAX (STRLEN/2)
2813 int nelem, i, j, k, nr;
2814 char *names[EGP_MAX];
2818 gnames = groups->grpname;
2820 nelem = str_nelem(val, EGP_MAX, names);
2823 gmx_fatal(FARGS, "The number of groups for %s is odd", option);
2825 nr = groups->grps[egcENER].nr;
2827 for (i = 0; i < nelem/2; i++)
2831 gmx_strcasecmp(names[2*i], *(gnames[groups->grps[egcENER].nm_ind[j]])))
2837 gmx_fatal(FARGS, "%s in %s is not an energy group\n",
2838 names[2*i], option);
2842 gmx_strcasecmp(names[2*i+1], *(gnames[groups->grps[egcENER].nm_ind[k]])))
2848 gmx_fatal(FARGS, "%s in %s is not an energy group\n",
2849 names[2*i+1], option);
2851 if ((j < nr) && (k < nr))
2853 ir->opts.egp_flags[nr*j+k] |= flag;
2854 ir->opts.egp_flags[nr*k+j] |= flag;
2862 void do_index(const char* mdparin, const char *ndx,
2865 t_inputrec *ir, rvec *v,
2869 gmx_groups_t *groups;
2873 char warnbuf[STRLEN], **gnames;
2874 int nr, ntcg, ntau_t, nref_t, nacc, nofg, nSA, nSA_points, nSA_time, nSA_temp;
2877 int nacg, nfreeze, nfrdim, nenergy, nvcm, nuser;
2878 char *ptr1[MAXPTR], *ptr2[MAXPTR], *ptr3[MAXPTR];
2879 int i, j, k, restnm;
2881 gmx_bool bExcl, bTable, bSetTCpar, bAnneal, bRest;
2882 int nQMmethod, nQMbasis, nQMcharge, nQMmult, nbSH, nCASorb, nCASelec,
2883 nSAon, nSAoff, nSAsteps, nQMg, nbOPT, nbTS;
2884 char warn_buf[STRLEN];
2888 fprintf(stderr, "processing index file...\n");
2894 snew(grps->index, 1);
2896 atoms_all = gmx_mtop_global_atoms(mtop);
2897 analyse(&atoms_all, grps, &gnames, FALSE, TRUE);
2898 free_t_atoms(&atoms_all, FALSE);
2902 grps = init_index(ndx, &gnames);
2905 groups = &mtop->groups;
2906 natoms = mtop->natoms;
2907 symtab = &mtop->symtab;
2909 snew(groups->grpname, grps->nr+1);
2911 for (i = 0; (i < grps->nr); i++)
2913 groups->grpname[i] = put_symtab(symtab, gnames[i]);
2915 groups->grpname[i] = put_symtab(symtab, "rest");
2917 srenew(gnames, grps->nr+1);
2918 gnames[restnm] = *(groups->grpname[i]);
2919 groups->ngrpname = grps->nr+1;
2921 set_warning_line(wi, mdparin, -1);
2923 ntau_t = str_nelem(tau_t, MAXPTR, ptr1);
2924 nref_t = str_nelem(ref_t, MAXPTR, ptr2);
2925 ntcg = str_nelem(tcgrps, MAXPTR, ptr3);
2926 if ((ntau_t != ntcg) || (nref_t != ntcg))
2928 gmx_fatal(FARGS, "Invalid T coupling input: %d groups, %d ref-t values and "
2929 "%d tau-t values", ntcg, nref_t, ntau_t);
2932 bSetTCpar = (ir->etc || EI_SD(ir->eI) || ir->eI == eiBD || EI_TPI(ir->eI));
2933 do_numbering(natoms, groups, ntcg, ptr3, grps, gnames, egcTC,
2934 restnm, bSetTCpar ? egrptpALL : egrptpALL_GENREST, bVerbose, wi);
2935 nr = groups->grps[egcTC].nr;
2937 snew(ir->opts.nrdf, nr);
2938 snew(ir->opts.tau_t, nr);
2939 snew(ir->opts.ref_t, nr);
2940 if (ir->eI == eiBD && ir->bd_fric == 0)
2942 fprintf(stderr, "bd-fric=0, so tau-t will be used as the inverse friction constant(s)\n");
2949 gmx_fatal(FARGS, "Not enough ref-t and tau-t values!");
2953 for (i = 0; (i < nr); i++)
2955 ir->opts.tau_t[i] = strtod(ptr1[i], NULL);
2956 if ((ir->eI == eiBD || ir->eI == eiSD2) && ir->opts.tau_t[i] <= 0)
2958 sprintf(warn_buf, "With integrator %s tau-t should be larger than 0", ei_names[ir->eI]);
2959 warning_error(wi, warn_buf);
2962 if (ir->etc != etcVRESCALE && ir->opts.tau_t[i] == 0)
2964 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.");
2967 if (ir->opts.tau_t[i] >= 0)
2969 tau_min = min(tau_min, ir->opts.tau_t[i]);
2972 if (ir->etc != etcNO && ir->nsttcouple == -1)
2974 ir->nsttcouple = ir_optimal_nsttcouple(ir);
2979 if ((ir->etc == etcNOSEHOOVER) && (ir->epc == epcBERENDSEN))
2981 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");
2983 if ((ir->epc == epcMTTK) && (ir->etc > etcNO))
2985 if (ir->nstpcouple != ir->nsttcouple)
2987 int mincouple = min(ir->nstpcouple, ir->nsttcouple);
2988 ir->nstpcouple = ir->nsttcouple = mincouple;
2989 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);
2990 warning_note(wi, warn_buf);
2994 /* velocity verlet with averaged kinetic energy KE = 0.5*(v(t+1/2) - v(t-1/2)) is implemented
2995 primarily for testing purposes, and does not work with temperature coupling other than 1 */
2997 if (ETC_ANDERSEN(ir->etc))
2999 if (ir->nsttcouple != 1)
3002 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");
3003 warning_note(wi, warn_buf);
3006 nstcmin = tcouple_min_integration_steps(ir->etc);
3009 if (tau_min/(ir->delta_t*ir->nsttcouple) < nstcmin)
3011 sprintf(warn_buf, "For proper integration of the %s thermostat, tau-t (%g) should be at least %d times larger than nsttcouple*dt (%g)",
3012 ETCOUPLTYPE(ir->etc),
3014 ir->nsttcouple*ir->delta_t);
3015 warning(wi, warn_buf);
3018 for (i = 0; (i < nr); i++)
3020 ir->opts.ref_t[i] = strtod(ptr2[i], NULL);
3021 if (ir->opts.ref_t[i] < 0)
3023 gmx_fatal(FARGS, "ref-t for group %d negative", i);
3026 /* set the lambda mc temperature to the md integrator temperature (which should be defined
3027 if we are in this conditional) if mc_temp is negative */
3028 if (ir->expandedvals->mc_temp < 0)
3030 ir->expandedvals->mc_temp = ir->opts.ref_t[0]; /*for now, set to the first reft */
3034 /* Simulated annealing for each group. There are nr groups */
3035 nSA = str_nelem(anneal, MAXPTR, ptr1);
3036 if (nSA == 1 && (ptr1[0][0] == 'n' || ptr1[0][0] == 'N'))
3040 if (nSA > 0 && nSA != nr)
3042 gmx_fatal(FARGS, "Not enough annealing values: %d (for %d groups)\n", nSA, nr);
3046 snew(ir->opts.annealing, nr);
3047 snew(ir->opts.anneal_npoints, nr);
3048 snew(ir->opts.anneal_time, nr);
3049 snew(ir->opts.anneal_temp, nr);
3050 for (i = 0; i < nr; i++)
3052 ir->opts.annealing[i] = eannNO;
3053 ir->opts.anneal_npoints[i] = 0;
3054 ir->opts.anneal_time[i] = NULL;
3055 ir->opts.anneal_temp[i] = NULL;
3060 for (i = 0; i < nr; i++)
3062 if (ptr1[i][0] == 'n' || ptr1[i][0] == 'N')
3064 ir->opts.annealing[i] = eannNO;
3066 else if (ptr1[i][0] == 's' || ptr1[i][0] == 'S')
3068 ir->opts.annealing[i] = eannSINGLE;
3071 else if (ptr1[i][0] == 'p' || ptr1[i][0] == 'P')
3073 ir->opts.annealing[i] = eannPERIODIC;
3079 /* Read the other fields too */
3080 nSA_points = str_nelem(anneal_npoints, MAXPTR, ptr1);
3081 if (nSA_points != nSA)
3083 gmx_fatal(FARGS, "Found %d annealing-npoints values for %d groups\n", nSA_points, nSA);
3085 for (k = 0, i = 0; i < nr; i++)
3087 ir->opts.anneal_npoints[i] = strtol(ptr1[i], NULL, 10);
3088 if (ir->opts.anneal_npoints[i] == 1)
3090 gmx_fatal(FARGS, "Please specify at least a start and an end point for annealing\n");
3092 snew(ir->opts.anneal_time[i], ir->opts.anneal_npoints[i]);
3093 snew(ir->opts.anneal_temp[i], ir->opts.anneal_npoints[i]);
3094 k += ir->opts.anneal_npoints[i];
3097 nSA_time = str_nelem(anneal_time, MAXPTR, ptr1);
3100 gmx_fatal(FARGS, "Found %d annealing-time values, wanter %d\n", nSA_time, k);
3102 nSA_temp = str_nelem(anneal_temp, MAXPTR, ptr2);
3105 gmx_fatal(FARGS, "Found %d annealing-temp values, wanted %d\n", nSA_temp, k);
3108 for (i = 0, k = 0; i < nr; i++)
3111 for (j = 0; j < ir->opts.anneal_npoints[i]; j++)
3113 ir->opts.anneal_time[i][j] = strtod(ptr1[k], NULL);
3114 ir->opts.anneal_temp[i][j] = strtod(ptr2[k], NULL);
3117 if (ir->opts.anneal_time[i][0] > (ir->init_t+GMX_REAL_EPS))
3119 gmx_fatal(FARGS, "First time point for annealing > init_t.\n");
3125 if (ir->opts.anneal_time[i][j] < ir->opts.anneal_time[i][j-1])
3127 gmx_fatal(FARGS, "Annealing timepoints out of order: t=%f comes after t=%f\n",
3128 ir->opts.anneal_time[i][j], ir->opts.anneal_time[i][j-1]);
3131 if (ir->opts.anneal_temp[i][j] < 0)
3133 gmx_fatal(FARGS, "Found negative temperature in annealing: %f\n", ir->opts.anneal_temp[i][j]);
3138 /* Print out some summary information, to make sure we got it right */
3139 for (i = 0, k = 0; i < nr; i++)
3141 if (ir->opts.annealing[i] != eannNO)
3143 j = groups->grps[egcTC].nm_ind[i];
3144 fprintf(stderr, "Simulated annealing for group %s: %s, %d timepoints\n",
3145 *(groups->grpname[j]), eann_names[ir->opts.annealing[i]],
3146 ir->opts.anneal_npoints[i]);
3147 fprintf(stderr, "Time (ps) Temperature (K)\n");
3148 /* All terms except the last one */
3149 for (j = 0; j < (ir->opts.anneal_npoints[i]-1); j++)
3151 fprintf(stderr, "%9.1f %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3154 /* Finally the last one */
3155 j = ir->opts.anneal_npoints[i]-1;
3156 if (ir->opts.annealing[i] == eannSINGLE)
3158 fprintf(stderr, "%9.1f- %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3162 fprintf(stderr, "%9.1f %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3163 if (fabs(ir->opts.anneal_temp[i][j]-ir->opts.anneal_temp[i][0]) > GMX_REAL_EPS)
3165 warning_note(wi, "There is a temperature jump when your annealing loops back.\n");
3174 if (ir->ePull != epullNO)
3176 make_pull_groups(ir->pull, pull_grp, grps, gnames);
3181 make_rotation_groups(ir->rot, rot_grp, grps, gnames);
3184 nacc = str_nelem(acc, MAXPTR, ptr1);
3185 nacg = str_nelem(accgrps, MAXPTR, ptr2);
3186 if (nacg*DIM != nacc)
3188 gmx_fatal(FARGS, "Invalid Acceleration input: %d groups and %d acc. values",
3191 do_numbering(natoms, groups, nacg, ptr2, grps, gnames, egcACC,
3192 restnm, egrptpALL_GENREST, bVerbose, wi);
3193 nr = groups->grps[egcACC].nr;
3194 snew(ir->opts.acc, nr);
3195 ir->opts.ngacc = nr;
3197 for (i = k = 0; (i < nacg); i++)
3199 for (j = 0; (j < DIM); j++, k++)
3201 ir->opts.acc[i][j] = strtod(ptr1[k], NULL);
3204 for (; (i < nr); i++)
3206 for (j = 0; (j < DIM); j++)
3208 ir->opts.acc[i][j] = 0;
3212 nfrdim = str_nelem(frdim, MAXPTR, ptr1);
3213 nfreeze = str_nelem(freeze, MAXPTR, ptr2);
3214 if (nfrdim != DIM*nfreeze)
3216 gmx_fatal(FARGS, "Invalid Freezing input: %d groups and %d freeze values",
3219 do_numbering(natoms, groups, nfreeze, ptr2, grps, gnames, egcFREEZE,
3220 restnm, egrptpALL_GENREST, bVerbose, wi);
3221 nr = groups->grps[egcFREEZE].nr;
3222 ir->opts.ngfrz = nr;
3223 snew(ir->opts.nFreeze, nr);
3224 for (i = k = 0; (i < nfreeze); i++)
3226 for (j = 0; (j < DIM); j++, k++)
3228 ir->opts.nFreeze[i][j] = (gmx_strncasecmp(ptr1[k], "Y", 1) == 0);
3229 if (!ir->opts.nFreeze[i][j])
3231 if (gmx_strncasecmp(ptr1[k], "N", 1) != 0)
3233 sprintf(warnbuf, "Please use Y(ES) or N(O) for freezedim only "
3234 "(not %s)", ptr1[k]);
3235 warning(wi, warn_buf);
3240 for (; (i < nr); i++)
3242 for (j = 0; (j < DIM); j++)
3244 ir->opts.nFreeze[i][j] = 0;
3248 nenergy = str_nelem(energy, MAXPTR, ptr1);
3249 do_numbering(natoms, groups, nenergy, ptr1, grps, gnames, egcENER,
3250 restnm, egrptpALL_GENREST, bVerbose, wi);
3251 add_wall_energrps(groups, ir->nwall, symtab);
3252 ir->opts.ngener = groups->grps[egcENER].nr;
3253 nvcm = str_nelem(vcm, MAXPTR, ptr1);
3255 do_numbering(natoms, groups, nvcm, ptr1, grps, gnames, egcVCM,
3256 restnm, nvcm == 0 ? egrptpALL_GENREST : egrptpPART, bVerbose, wi);
3259 warning(wi, "Some atoms are not part of any center of mass motion removal group.\n"
3260 "This may lead to artifacts.\n"
3261 "In most cases one should use one group for the whole system.");
3264 /* Now we have filled the freeze struct, so we can calculate NRDF */
3265 calc_nrdf(mtop, ir, gnames);
3271 /* Must check per group! */
3272 for (i = 0; (i < ir->opts.ngtc); i++)
3274 ntot += ir->opts.nrdf[i];
3276 if (ntot != (DIM*natoms))
3278 fac = sqrt(ntot/(DIM*natoms));
3281 fprintf(stderr, "Scaling velocities by a factor of %.3f to account for constraints\n"
3282 "and removal of center of mass motion\n", fac);
3284 for (i = 0; (i < natoms); i++)
3286 svmul(fac, v[i], v[i]);
3291 nuser = str_nelem(user1, MAXPTR, ptr1);
3292 do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcUser1,
3293 restnm, egrptpALL_GENREST, bVerbose, wi);
3294 nuser = str_nelem(user2, MAXPTR, ptr1);
3295 do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcUser2,
3296 restnm, egrptpALL_GENREST, bVerbose, wi);
3297 nuser = str_nelem(xtc_grps, MAXPTR, ptr1);
3298 do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcXTC,
3299 restnm, egrptpONE, bVerbose, wi);
3300 nofg = str_nelem(orirefitgrp, MAXPTR, ptr1);
3301 do_numbering(natoms, groups, nofg, ptr1, grps, gnames, egcORFIT,
3302 restnm, egrptpALL_GENREST, bVerbose, wi);
3304 /* QMMM input processing */
3305 nQMg = str_nelem(QMMM, MAXPTR, ptr1);
3306 nQMmethod = str_nelem(QMmethod, MAXPTR, ptr2);
3307 nQMbasis = str_nelem(QMbasis, MAXPTR, ptr3);
3308 if ((nQMmethod != nQMg) || (nQMbasis != nQMg))
3310 gmx_fatal(FARGS, "Invalid QMMM input: %d groups %d basissets"
3311 " and %d methods\n", nQMg, nQMbasis, nQMmethod);
3313 /* group rest, if any, is always MM! */
3314 do_numbering(natoms, groups, nQMg, ptr1, grps, gnames, egcQMMM,
3315 restnm, egrptpALL_GENREST, bVerbose, wi);
3316 nr = nQMg; /*atoms->grps[egcQMMM].nr;*/
3317 ir->opts.ngQM = nQMg;
3318 snew(ir->opts.QMmethod, nr);
3319 snew(ir->opts.QMbasis, nr);
3320 for (i = 0; i < nr; i++)
3322 /* input consists of strings: RHF CASSCF PM3 .. These need to be
3323 * converted to the corresponding enum in names.c
3325 ir->opts.QMmethod[i] = search_QMstring(ptr2[i], eQMmethodNR,
3327 ir->opts.QMbasis[i] = search_QMstring(ptr3[i], eQMbasisNR,
3331 nQMmult = str_nelem(QMmult, MAXPTR, ptr1);
3332 nQMcharge = str_nelem(QMcharge, MAXPTR, ptr2);
3333 nbSH = str_nelem(bSH, MAXPTR, ptr3);
3334 snew(ir->opts.QMmult, nr);
3335 snew(ir->opts.QMcharge, nr);
3336 snew(ir->opts.bSH, nr);
3338 for (i = 0; i < nr; i++)
3340 ir->opts.QMmult[i] = strtol(ptr1[i], NULL, 10);
3341 ir->opts.QMcharge[i] = strtol(ptr2[i], NULL, 10);
3342 ir->opts.bSH[i] = (gmx_strncasecmp(ptr3[i], "Y", 1) == 0);
3345 nCASelec = str_nelem(CASelectrons, MAXPTR, ptr1);
3346 nCASorb = str_nelem(CASorbitals, MAXPTR, ptr2);
3347 snew(ir->opts.CASelectrons, nr);
3348 snew(ir->opts.CASorbitals, nr);
3349 for (i = 0; i < nr; i++)
3351 ir->opts.CASelectrons[i] = strtol(ptr1[i], NULL, 10);
3352 ir->opts.CASorbitals[i] = strtol(ptr2[i], NULL, 10);
3354 /* special optimization options */
3356 nbOPT = str_nelem(bOPT, MAXPTR, ptr1);
3357 nbTS = str_nelem(bTS, MAXPTR, ptr2);
3358 snew(ir->opts.bOPT, nr);
3359 snew(ir->opts.bTS, nr);
3360 for (i = 0; i < nr; i++)
3362 ir->opts.bOPT[i] = (gmx_strncasecmp(ptr1[i], "Y", 1) == 0);
3363 ir->opts.bTS[i] = (gmx_strncasecmp(ptr2[i], "Y", 1) == 0);
3365 nSAon = str_nelem(SAon, MAXPTR, ptr1);
3366 nSAoff = str_nelem(SAoff, MAXPTR, ptr2);
3367 nSAsteps = str_nelem(SAsteps, MAXPTR, ptr3);
3368 snew(ir->opts.SAon, nr);
3369 snew(ir->opts.SAoff, nr);
3370 snew(ir->opts.SAsteps, nr);
3372 for (i = 0; i < nr; i++)
3374 ir->opts.SAon[i] = strtod(ptr1[i], NULL);
3375 ir->opts.SAoff[i] = strtod(ptr2[i], NULL);
3376 ir->opts.SAsteps[i] = strtol(ptr3[i], NULL, 10);
3378 /* end of QMMM input */
3382 for (i = 0; (i < egcNR); i++)
3384 fprintf(stderr, "%-16s has %d element(s):", gtypes[i], groups->grps[i].nr);
3385 for (j = 0; (j < groups->grps[i].nr); j++)
3387 fprintf(stderr, " %s", *(groups->grpname[groups->grps[i].nm_ind[j]]));
3389 fprintf(stderr, "\n");
3393 nr = groups->grps[egcENER].nr;
3394 snew(ir->opts.egp_flags, nr*nr);
3396 bExcl = do_egp_flag(ir, groups, "energygrp-excl", egpexcl, EGP_EXCL);
3397 if (bExcl && ir->cutoff_scheme == ecutsVERLET)
3399 warning_error(wi, "Energy group exclusions are not (yet) implemented for the Verlet scheme");
3401 if (bExcl && EEL_FULL(ir->coulombtype))
3403 warning(wi, "Can not exclude the lattice Coulomb energy between energy groups");
3406 bTable = do_egp_flag(ir, groups, "energygrp-table", egptable, EGP_TABLE);
3407 if (bTable && !(ir->vdwtype == evdwUSER) &&
3408 !(ir->coulombtype == eelUSER) && !(ir->coulombtype == eelPMEUSER) &&
3409 !(ir->coulombtype == eelPMEUSERSWITCH))
3411 gmx_fatal(FARGS, "Can only have energy group pair tables in combination with user tables for VdW and/or Coulomb");
3414 decode_cos(efield_x, &(ir->ex[XX]));
3415 decode_cos(efield_xt, &(ir->et[XX]));
3416 decode_cos(efield_y, &(ir->ex[YY]));
3417 decode_cos(efield_yt, &(ir->et[YY]));
3418 decode_cos(efield_z, &(ir->ex[ZZ]));
3419 decode_cos(efield_zt, &(ir->et[ZZ]));
3423 do_adress_index(ir->adress, groups, gnames, &(ir->opts), wi);
3426 for (i = 0; (i < grps->nr); i++)
3438 static void check_disre(gmx_mtop_t *mtop)
3440 gmx_ffparams_t *ffparams;
3441 t_functype *functype;
3443 int i, ndouble, ftype;
3444 int label, old_label;
3446 if (gmx_mtop_ftype_count(mtop, F_DISRES) > 0)
3448 ffparams = &mtop->ffparams;
3449 functype = ffparams->functype;
3450 ip = ffparams->iparams;
3453 for (i = 0; i < ffparams->ntypes; i++)
3455 ftype = functype[i];
3456 if (ftype == F_DISRES)
3458 label = ip[i].disres.label;
3459 if (label == old_label)
3461 fprintf(stderr, "Distance restraint index %d occurs twice\n", label);
3469 gmx_fatal(FARGS, "Found %d double distance restraint indices,\n"
3470 "probably the parameters for multiple pairs in one restraint "
3471 "are not identical\n", ndouble);
3476 static gmx_bool absolute_reference(t_inputrec *ir, gmx_mtop_t *sys,
3477 gmx_bool posres_only,
3481 gmx_mtop_ilistloop_t iloop;
3491 for (d = 0; d < DIM; d++)
3493 AbsRef[d] = (d < ndof_com(ir) ? 0 : 1);
3495 /* Check for freeze groups */
3496 for (g = 0; g < ir->opts.ngfrz; g++)
3498 for (d = 0; d < DIM; d++)
3500 if (ir->opts.nFreeze[g][d] != 0)
3508 /* Check for position restraints */
3509 iloop = gmx_mtop_ilistloop_init(sys);
3510 while (gmx_mtop_ilistloop_next(iloop, &ilist, &nmol))
3513 (AbsRef[XX] == 0 || AbsRef[YY] == 0 || AbsRef[ZZ] == 0))
3515 for (i = 0; i < ilist[F_POSRES].nr; i += 2)
3517 pr = &sys->ffparams.iparams[ilist[F_POSRES].iatoms[i]];
3518 for (d = 0; d < DIM; d++)
3520 if (pr->posres.fcA[d] != 0)
3526 for (i = 0; i < ilist[F_FBPOSRES].nr; i += 2)
3528 /* Check for flat-bottom posres */
3529 pr = &sys->ffparams.iparams[ilist[F_FBPOSRES].iatoms[i]];
3530 if (pr->fbposres.k != 0)
3532 switch (pr->fbposres.geom)
3534 case efbposresSPHERE:
3535 AbsRef[XX] = AbsRef[YY] = AbsRef[ZZ] = 1;
3537 case efbposresCYLINDER:
3538 AbsRef[XX] = AbsRef[YY] = 1;
3540 case efbposresX: /* d=XX */
3541 case efbposresY: /* d=YY */
3542 case efbposresZ: /* d=ZZ */
3543 d = pr->fbposres.geom - efbposresX;
3547 gmx_fatal(FARGS, " Invalid geometry for flat-bottom position restraint.\n"
3548 "Expected nr between 1 and %d. Found %d\n", efbposresNR-1,
3556 return (AbsRef[XX] != 0 && AbsRef[YY] != 0 && AbsRef[ZZ] != 0);
3560 check_combination_rule_differences(const gmx_mtop_t *mtop, int state,
3561 gmx_bool *bC6ParametersWorkWithGeometricRules,
3562 gmx_bool *bC6ParametersWorkWithLBRules,
3563 gmx_bool *bLBRulesPossible)
3565 int ntypes, tpi, tpj, thisLBdiff, thisgeomdiff;
3568 #if (defined SIZEOF_LONG_LONG_INT) && (SIZEOF_LONG_LONG_INT >= 8)
3569 long long int npair, npair_ij;
3571 double npair, npair_ij;
3573 double geometricdiff, LBdiff;
3574 double c6i, c6j, c12i, c12j;
3575 double c6, c6_geometric, c6_LB;
3576 double sigmai, sigmaj, epsi, epsj;
3577 gmx_bool bCanDoLBRules, bCanDoGeometricRules;
3580 /* A tolerance of 1e-5 seems reasonable for (possibly hand-typed)
3581 * force-field floating point parameters.
3584 ptr = getenv("GMX_LJCOMB_TOL");
3589 sscanf(ptr, "%lf", &dbl);
3593 *bC6ParametersWorkWithLBRules = TRUE;
3594 *bC6ParametersWorkWithGeometricRules = TRUE;
3595 bCanDoLBRules = TRUE;
3596 bCanDoGeometricRules = TRUE;
3598 ntypes = mtop->ffparams.atnr;
3599 snew(typecount, ntypes);
3600 gmx_mtop_count_atomtypes(mtop, state, typecount);
3601 geometricdiff = LBdiff = 0.0;
3602 *bLBRulesPossible = TRUE;
3603 for (tpi = 0; tpi < ntypes; ++tpi)
3605 c6i = mtop->ffparams.iparams[(ntypes + 1) * tpi].lj.c6;
3606 c12i = mtop->ffparams.iparams[(ntypes + 1) * tpi].lj.c12;
3607 for (tpj = tpi; tpj < ntypes; ++tpj)
3609 c6j = mtop->ffparams.iparams[(ntypes + 1) * tpj].lj.c6;
3610 c12j = mtop->ffparams.iparams[(ntypes + 1) * tpj].lj.c12;
3611 c6 = mtop->ffparams.iparams[ntypes * tpi + tpj].lj.c6;
3612 c6_geometric = sqrt(c6i * c6j);
3613 if (!gmx_numzero(c6_geometric))
3615 if (!gmx_numzero(c12i) && !gmx_numzero(c12j))
3617 sigmai = pow(c12i / c6i, 1.0/6.0);
3618 sigmaj = pow(c12j / c6j, 1.0/6.0);
3619 epsi = c6i * c6i /(4.0 * c12i);
3620 epsj = c6j * c6j /(4.0 * c12j);
3621 c6_LB = 4.0 * pow(epsi * epsj, 1.0/2.0) * pow(0.5 * (sigmai + sigmaj), 6);
3625 *bLBRulesPossible = FALSE;
3626 c6_LB = c6_geometric;
3628 bCanDoLBRules = gmx_within_tol(c6_LB, c6, tol);
3631 if (FALSE == bCanDoLBRules)
3633 *bC6ParametersWorkWithLBRules = FALSE;
3636 bCanDoGeometricRules = gmx_within_tol(c6_geometric, c6, tol);
3638 if (FALSE == bCanDoGeometricRules)
3640 *bC6ParametersWorkWithGeometricRules = FALSE;
3648 check_combination_rules(const t_inputrec *ir, const gmx_mtop_t *mtop,
3652 gmx_bool bLBRulesPossible, bC6ParametersWorkWithGeometricRules, bC6ParametersWorkWithLBRules;
3654 check_combination_rule_differences(mtop, 0,
3655 &bC6ParametersWorkWithGeometricRules,
3656 &bC6ParametersWorkWithLBRules,
3658 if (ir->ljpme_combination_rule == eljpmeLB)
3660 if (FALSE == bC6ParametersWorkWithLBRules || FALSE == bLBRulesPossible)
3662 warning(wi, "You are using arithmetic-geometric combination rules "
3663 "in LJ-PME, but your non-bonded C6 parameters do not "
3664 "follow these rules.");
3669 if (FALSE == bC6ParametersWorkWithGeometricRules)
3671 if (ir->eDispCorr != edispcNO)
3673 warning_note(wi, "You are using geometric combination rules in "
3674 "LJ-PME, but your non-bonded C6 parameters do "
3675 "not follow these rules. "
3676 "If your force field uses Lorentz-Berthelot combination rules this "
3677 "will introduce small errors in the forces and energies in "
3678 "your simulations. Dispersion correction will correct total "
3679 "energy and/or pressure, but not forces or surface tensions. "
3680 "Please check the LJ-PME section in the manual "
3681 "before proceeding further.");
3685 warning_note(wi, "You are using geometric combination rules in "
3686 "LJ-PME, but your non-bonded C6 parameters do "
3687 "not follow these rules. "
3688 "If your force field uses Lorentz-Berthelot combination rules this "
3689 "will introduce small errors in the forces and energies in "
3690 "your simulations. Consider using dispersion correction "
3691 "for the total energy and pressure. "
3692 "Please check the LJ-PME section in the manual "
3693 "before proceeding further.");
3699 void triple_check(const char *mdparin, t_inputrec *ir, gmx_mtop_t *sys,
3703 int i, m, g, nmol, npct;
3704 gmx_bool bCharge, bAcc;
3705 real gdt_max, *mgrp, mt;
3707 gmx_mtop_atomloop_block_t aloopb;
3708 gmx_mtop_atomloop_all_t aloop;
3711 char warn_buf[STRLEN];
3713 set_warning_line(wi, mdparin, -1);
3715 if (EI_DYNAMICS(ir->eI) && !EI_SD(ir->eI) && ir->eI != eiBD &&
3716 ir->comm_mode == ecmNO &&
3717 !(absolute_reference(ir, sys, FALSE, AbsRef) || ir->nsteps <= 10))
3719 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");
3722 /* Check for pressure coupling with absolute position restraints */
3723 if (ir->epc != epcNO && ir->refcoord_scaling == erscNO)
3725 absolute_reference(ir, sys, TRUE, AbsRef);
3727 for (m = 0; m < DIM; m++)
3729 if (AbsRef[m] && norm2(ir->compress[m]) > 0)
3731 warning(wi, "You are using pressure coupling with absolute position restraints, this will give artifacts. Use the refcoord_scaling option.");
3739 aloopb = gmx_mtop_atomloop_block_init(sys);
3740 while (gmx_mtop_atomloop_block_next(aloopb, &atom, &nmol))
3742 if (atom->q != 0 || atom->qB != 0)
3750 if (EEL_FULL(ir->coulombtype))
3753 "You are using full electrostatics treatment %s for a system without charges.\n"
3754 "This costs a lot of performance for just processing zeros, consider using %s instead.\n",
3755 EELTYPE(ir->coulombtype), EELTYPE(eelCUT));
3756 warning(wi, err_buf);
3761 if (ir->coulombtype == eelCUT && ir->rcoulomb > 0 && !ir->implicit_solvent)
3764 "You are using a plain Coulomb cut-off, which might produce artifacts.\n"
3765 "You might want to consider using %s electrostatics.\n",
3767 warning_note(wi, err_buf);
3771 /* Check if combination rules used in LJ-PME are the same as in the force field */
3772 if (EVDW_PME(ir->vdwtype))
3774 check_combination_rules(ir, sys, wi);
3777 /* Generalized reaction field */
3778 if (ir->opts.ngtc == 0)
3780 sprintf(err_buf, "No temperature coupling while using coulombtype %s",
3782 CHECK(ir->coulombtype == eelGRF);
3786 sprintf(err_buf, "When using coulombtype = %s"
3787 " ref-t for temperature coupling should be > 0",
3789 CHECK((ir->coulombtype == eelGRF) && (ir->opts.ref_t[0] <= 0));
3792 if (ir->eI == eiSD1 &&
3793 (gmx_mtop_ftype_count(sys, F_CONSTR) > 0 ||
3794 gmx_mtop_ftype_count(sys, F_SETTLE) > 0))
3796 sprintf(warn_buf, "With constraints integrator %s is less accurate, consider using %s instead", ei_names[ir->eI], ei_names[eiSD2]);
3797 warning_note(wi, warn_buf);
3801 for (i = 0; (i < sys->groups.grps[egcACC].nr); i++)
3803 for (m = 0; (m < DIM); m++)
3805 if (fabs(ir->opts.acc[i][m]) > 1e-6)
3814 snew(mgrp, sys->groups.grps[egcACC].nr);
3815 aloop = gmx_mtop_atomloop_all_init(sys);
3816 while (gmx_mtop_atomloop_all_next(aloop, &i, &atom))
3818 mgrp[ggrpnr(&sys->groups, egcACC, i)] += atom->m;
3821 for (i = 0; (i < sys->groups.grps[egcACC].nr); i++)
3823 for (m = 0; (m < DIM); m++)
3825 acc[m] += ir->opts.acc[i][m]*mgrp[i];
3829 for (m = 0; (m < DIM); m++)
3831 if (fabs(acc[m]) > 1e-6)
3833 const char *dim[DIM] = { "X", "Y", "Z" };
3835 "Net Acceleration in %s direction, will %s be corrected\n",
3836 dim[m], ir->nstcomm != 0 ? "" : "not");
3837 if (ir->nstcomm != 0 && m < ndof_com(ir))
3840 for (i = 0; (i < sys->groups.grps[egcACC].nr); i++)
3842 ir->opts.acc[i][m] -= acc[m];
3850 if (ir->efep != efepNO && ir->fepvals->sc_alpha != 0 &&
3851 !gmx_within_tol(sys->ffparams.reppow, 12.0, 10*GMX_DOUBLE_EPS))
3853 gmx_fatal(FARGS, "Soft-core interactions are only supported with VdW repulsion power 12");
3856 if (ir->ePull != epullNO)
3858 if (ir->pull->grp[0].nat == 0)
3860 absolute_reference(ir, sys, FALSE, AbsRef);
3861 for (m = 0; m < DIM; m++)
3863 if (ir->pull->dim[m] && !AbsRef[m])
3865 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.");
3871 if (ir->pull->eGeom == epullgDIRPBC)
3873 for (i = 0; i < 3; i++)
3875 for (m = 0; m <= i; m++)
3877 if ((ir->epc != epcNO && ir->compress[i][m] != 0) ||
3878 ir->deform[i][m] != 0)
3880 for (g = 1; g < ir->pull->ngrp; g++)
3882 if (ir->pull->grp[g].vec[m] != 0)
3884 gmx_fatal(FARGS, "Can not have dynamic box while using pull geometry '%s' (dim %c)", EPULLGEOM(ir->pull->eGeom), 'x'+m);
3896 void double_check(t_inputrec *ir, matrix box, gmx_bool bConstr, warninp_t wi)
3900 char warn_buf[STRLEN];
3903 ptr = check_box(ir->ePBC, box);
3906 warning_error(wi, ptr);
3909 if (bConstr && ir->eConstrAlg == econtSHAKE)
3911 if (ir->shake_tol <= 0.0)
3913 sprintf(warn_buf, "ERROR: shake-tol must be > 0 instead of %g\n",
3915 warning_error(wi, warn_buf);
3918 if (IR_TWINRANGE(*ir) && ir->nstlist > 1)
3920 sprintf(warn_buf, "With twin-range cut-off's and SHAKE the virial and the pressure are incorrect.");
3921 if (ir->epc == epcNO)
3923 warning(wi, warn_buf);
3927 warning_error(wi, warn_buf);
3932 if ( (ir->eConstrAlg == econtLINCS) && bConstr)
3934 /* If we have Lincs constraints: */
3935 if (ir->eI == eiMD && ir->etc == etcNO &&
3936 ir->eConstrAlg == econtLINCS && ir->nLincsIter == 1)
3938 sprintf(warn_buf, "For energy conservation with LINCS, lincs_iter should be 2 or larger.\n");
3939 warning_note(wi, warn_buf);
3942 if ((ir->eI == eiCG || ir->eI == eiLBFGS) && (ir->nProjOrder < 8))
3944 sprintf(warn_buf, "For accurate %s with LINCS constraints, lincs-order should be 8 or more.", ei_names[ir->eI]);
3945 warning_note(wi, warn_buf);
3947 if (ir->epc == epcMTTK)
3949 warning_error(wi, "MTTK not compatible with lincs -- use shake instead.");
3953 if (ir->LincsWarnAngle > 90.0)
3955 sprintf(warn_buf, "lincs-warnangle can not be larger than 90 degrees, setting it to 90.\n");
3956 warning(wi, warn_buf);
3957 ir->LincsWarnAngle = 90.0;
3960 if (ir->ePBC != epbcNONE)
3962 if (ir->nstlist == 0)
3964 warning(wi, "With nstlist=0 atoms are only put into the box at step 0, therefore drifting atoms might cause the simulation to crash.");
3966 bTWIN = (ir->rlistlong > ir->rlist);
3967 if (ir->ns_type == ensGRID)
3969 if (sqr(ir->rlistlong) >= max_cutoff2(ir->ePBC, box))
3971 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",
3972 bTWIN ? (ir->rcoulomb == ir->rlistlong ? "rcoulomb" : "rvdw") : "rlist");
3973 warning_error(wi, warn_buf);
3978 min_size = min(box[XX][XX], min(box[YY][YY], box[ZZ][ZZ]));
3979 if (2*ir->rlistlong >= min_size)
3981 sprintf(warn_buf, "ERROR: One of the box lengths is smaller than twice the cut-off length. Increase the box size or decrease rlist.");
3982 warning_error(wi, warn_buf);
3985 fprintf(stderr, "Grid search might allow larger cut-off's than simple search with triclinic boxes.");
3992 void check_chargegroup_radii(const gmx_mtop_t *mtop, const t_inputrec *ir,
3996 real rvdw1, rvdw2, rcoul1, rcoul2;
3997 char warn_buf[STRLEN];
3999 calc_chargegroup_radii(mtop, x, &rvdw1, &rvdw2, &rcoul1, &rcoul2);
4003 printf("Largest charge group radii for Van der Waals: %5.3f, %5.3f nm\n",
4008 printf("Largest charge group radii for Coulomb: %5.3f, %5.3f nm\n",
4014 if (rvdw1 + rvdw2 > ir->rlist ||
4015 rcoul1 + rcoul2 > ir->rlist)
4018 "The sum of the two largest charge group radii (%f) "
4019 "is larger than rlist (%f)\n",
4020 max(rvdw1+rvdw2, rcoul1+rcoul2), ir->rlist);
4021 warning(wi, warn_buf);
4025 /* Here we do not use the zero at cut-off macro,
4026 * since user defined interactions might purposely
4027 * not be zero at the cut-off.
4029 if ((EVDW_IS_ZERO_AT_CUTOFF(ir->vdwtype) ||
4030 ir->vdw_modifier != eintmodNONE) &&
4031 rvdw1 + rvdw2 > ir->rlistlong - ir->rvdw)
4033 sprintf(warn_buf, "The sum of the two largest charge group "
4034 "radii (%f) is larger than %s (%f) - rvdw (%f).\n"
4035 "With exact cut-offs, better performance can be "
4036 "obtained with cutoff-scheme = %s, because it "
4037 "does not use charge groups at all.",
4039 ir->rlistlong > ir->rlist ? "rlistlong" : "rlist",
4040 ir->rlistlong, ir->rvdw,
4041 ecutscheme_names[ecutsVERLET]);
4044 warning(wi, warn_buf);
4048 warning_note(wi, warn_buf);
4051 if ((EEL_IS_ZERO_AT_CUTOFF(ir->coulombtype) ||
4052 ir->coulomb_modifier != eintmodNONE) &&
4053 rcoul1 + rcoul2 > ir->rlistlong - ir->rcoulomb)
4055 sprintf(warn_buf, "The sum of the two largest charge group radii (%f) is larger than %s (%f) - rcoulomb (%f).\n"
4056 "With exact cut-offs, better performance can be obtained with cutoff-scheme = %s, because it does not use charge groups at all.",
4058 ir->rlistlong > ir->rlist ? "rlistlong" : "rlist",
4059 ir->rlistlong, ir->rcoulomb,
4060 ecutscheme_names[ecutsVERLET]);
4063 warning(wi, warn_buf);
4067 warning_note(wi, warn_buf);