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_drift > 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_drift <= 0)
330 if (ir->verletbuf_drift == 0)
332 warning_error(wi, "Can not have an energy drift 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-drift > 0. Will set rlist using verlet-buffer-drift.");
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 drift with verlet-buffer-drift > 0. Either use temperature coupling or set rlist yourself together with verlet-buffer-drift = -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-drift > 0. You are using at least one unbounded dimension, so no volume can be computed. Either use a finite box, or set rlist yourself together with verlet-buffer-drift = -1.");
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");
1717 CCTYPE ("VARIOUS PREPROCESSING OPTIONS");
1718 CTYPE ("Preprocessor information: use cpp syntax.");
1719 CTYPE ("e.g.: -I/home/joe/doe -I/home/mary/roe");
1720 STYPE ("include", opts->include, NULL);
1721 CTYPE ("e.g.: -DPOSRES -DFLEXIBLE (note these variable names are case sensitive)");
1722 STYPE ("define", opts->define, NULL);
1724 CCTYPE ("RUN CONTROL PARAMETERS");
1725 EETYPE("integrator", ir->eI, ei_names);
1726 CTYPE ("Start time and timestep in ps");
1727 RTYPE ("tinit", ir->init_t, 0.0);
1728 RTYPE ("dt", ir->delta_t, 0.001);
1729 STEPTYPE ("nsteps", ir->nsteps, 0);
1730 CTYPE ("For exact run continuation or redoing part of a run");
1731 STEPTYPE ("init-step", ir->init_step, 0);
1732 CTYPE ("Part index is updated automatically on checkpointing (keeps files separate)");
1733 ITYPE ("simulation-part", ir->simulation_part, 1);
1734 CTYPE ("mode for center of mass motion removal");
1735 EETYPE("comm-mode", ir->comm_mode, ecm_names);
1736 CTYPE ("number of steps for center of mass motion removal");
1737 ITYPE ("nstcomm", ir->nstcomm, 100);
1738 CTYPE ("group(s) for center of mass motion removal");
1739 STYPE ("comm-grps", vcm, NULL);
1741 CCTYPE ("LANGEVIN DYNAMICS OPTIONS");
1742 CTYPE ("Friction coefficient (amu/ps) and random seed");
1743 RTYPE ("bd-fric", ir->bd_fric, 0.0);
1744 ITYPE ("ld-seed", ir->ld_seed, 1993);
1747 CCTYPE ("ENERGY MINIMIZATION OPTIONS");
1748 CTYPE ("Force tolerance and initial step-size");
1749 RTYPE ("emtol", ir->em_tol, 10.0);
1750 RTYPE ("emstep", ir->em_stepsize, 0.01);
1751 CTYPE ("Max number of iterations in relax-shells");
1752 ITYPE ("niter", ir->niter, 20);
1753 CTYPE ("Step size (ps^2) for minimization of flexible constraints");
1754 RTYPE ("fcstep", ir->fc_stepsize, 0);
1755 CTYPE ("Frequency of steepest descents steps when doing CG");
1756 ITYPE ("nstcgsteep", ir->nstcgsteep, 1000);
1757 ITYPE ("nbfgscorr", ir->nbfgscorr, 10);
1759 CCTYPE ("TEST PARTICLE INSERTION OPTIONS");
1760 RTYPE ("rtpi", ir->rtpi, 0.05);
1762 /* Output options */
1763 CCTYPE ("OUTPUT CONTROL OPTIONS");
1764 CTYPE ("Output frequency for coords (x), velocities (v) and forces (f)");
1765 ITYPE ("nstxout", ir->nstxout, 0);
1766 ITYPE ("nstvout", ir->nstvout, 0);
1767 ITYPE ("nstfout", ir->nstfout, 0);
1768 ir->nstcheckpoint = 1000;
1769 CTYPE ("Output frequency for energies to log file and energy file");
1770 ITYPE ("nstlog", ir->nstlog, 1000);
1771 ITYPE ("nstcalcenergy", ir->nstcalcenergy, 100);
1772 ITYPE ("nstenergy", ir->nstenergy, 1000);
1773 CTYPE ("Output frequency and precision for .xtc file");
1774 ITYPE ("nstxtcout", ir->nstxtcout, 0);
1775 RTYPE ("xtc-precision", ir->xtcprec, 1000.0);
1776 CTYPE ("This selects the subset of atoms for the .xtc file. You can");
1777 CTYPE ("select multiple groups. By default all atoms will be written.");
1778 STYPE ("xtc-grps", xtc_grps, NULL);
1779 CTYPE ("Selection of energy groups");
1780 STYPE ("energygrps", energy, NULL);
1782 /* Neighbor searching */
1783 CCTYPE ("NEIGHBORSEARCHING PARAMETERS");
1784 CTYPE ("cut-off scheme (group: using charge groups, Verlet: particle based cut-offs)");
1785 EETYPE("cutoff-scheme", ir->cutoff_scheme, ecutscheme_names);
1786 CTYPE ("nblist update frequency");
1787 ITYPE ("nstlist", ir->nstlist, 10);
1788 CTYPE ("ns algorithm (simple or grid)");
1789 EETYPE("ns-type", ir->ns_type, ens_names);
1790 /* set ndelta to the optimal value of 2 */
1792 CTYPE ("Periodic boundary conditions: xyz, no, xy");
1793 EETYPE("pbc", ir->ePBC, epbc_names);
1794 EETYPE("periodic-molecules", ir->bPeriodicMols, yesno_names);
1795 CTYPE ("Allowed energy drift due to the Verlet buffer in kJ/mol/ps per atom,");
1796 CTYPE ("a value of -1 means: use rlist");
1797 RTYPE("verlet-buffer-drift", ir->verletbuf_drift, 0.005);
1798 CTYPE ("nblist cut-off");
1799 RTYPE ("rlist", ir->rlist, 1.0);
1800 CTYPE ("long-range cut-off for switched potentials");
1801 RTYPE ("rlistlong", ir->rlistlong, -1);
1802 ITYPE ("nstcalclr", ir->nstcalclr, -1);
1804 /* Electrostatics */
1805 CCTYPE ("OPTIONS FOR ELECTROSTATICS AND VDW");
1806 CTYPE ("Method for doing electrostatics");
1807 EETYPE("coulombtype", ir->coulombtype, eel_names);
1808 EETYPE("coulomb-modifier", ir->coulomb_modifier, eintmod_names);
1809 CTYPE ("cut-off lengths");
1810 RTYPE ("rcoulomb-switch", ir->rcoulomb_switch, 0.0);
1811 RTYPE ("rcoulomb", ir->rcoulomb, 1.0);
1812 CTYPE ("Relative dielectric constant for the medium and the reaction field");
1813 RTYPE ("epsilon-r", ir->epsilon_r, 1.0);
1814 RTYPE ("epsilon-rf", ir->epsilon_rf, 0.0);
1815 CTYPE ("Method for doing Van der Waals");
1816 EETYPE("vdw-type", ir->vdwtype, evdw_names);
1817 EETYPE("vdw-modifier", ir->vdw_modifier, eintmod_names);
1818 CTYPE ("cut-off lengths");
1819 RTYPE ("rvdw-switch", ir->rvdw_switch, 0.0);
1820 RTYPE ("rvdw", ir->rvdw, 1.0);
1821 CTYPE ("Apply long range dispersion corrections for Energy and Pressure");
1822 EETYPE("DispCorr", ir->eDispCorr, edispc_names);
1823 CTYPE ("Extension of the potential lookup tables beyond the cut-off");
1824 RTYPE ("table-extension", ir->tabext, 1.0);
1825 CTYPE ("Separate tables between energy group pairs");
1826 STYPE ("energygrp-table", egptable, NULL);
1827 CTYPE ("Spacing for the PME/PPPM FFT grid");
1828 RTYPE ("fourierspacing", ir->fourier_spacing, 0.12);
1829 CTYPE ("FFT grid size, when a value is 0 fourierspacing will be used");
1830 ITYPE ("fourier-nx", ir->nkx, 0);
1831 ITYPE ("fourier-ny", ir->nky, 0);
1832 ITYPE ("fourier-nz", ir->nkz, 0);
1833 CTYPE ("EWALD/PME/PPPM parameters");
1834 ITYPE ("pme-order", ir->pme_order, 4);
1835 RTYPE ("ewald-rtol", ir->ewald_rtol, 0.00001);
1836 RTYPE ("ewald-rtol-lj", ir->ewald_rtol_lj, 0.001);
1837 EETYPE("lj-pme-comb-rule", ir->ljpme_combination_rule, eljpme_names);
1838 EETYPE("ewald-geometry", ir->ewald_geometry, eewg_names);
1839 RTYPE ("epsilon-surface", ir->epsilon_surface, 0.0);
1840 EETYPE("optimize-fft", ir->bOptFFT, yesno_names);
1842 CCTYPE("IMPLICIT SOLVENT ALGORITHM");
1843 EETYPE("implicit-solvent", ir->implicit_solvent, eis_names);
1845 CCTYPE ("GENERALIZED BORN ELECTROSTATICS");
1846 CTYPE ("Algorithm for calculating Born radii");
1847 EETYPE("gb-algorithm", ir->gb_algorithm, egb_names);
1848 CTYPE ("Frequency of calculating the Born radii inside rlist");
1849 ITYPE ("nstgbradii", ir->nstgbradii, 1);
1850 CTYPE ("Cutoff for Born radii calculation; the contribution from atoms");
1851 CTYPE ("between rlist and rgbradii is updated every nstlist steps");
1852 RTYPE ("rgbradii", ir->rgbradii, 1.0);
1853 CTYPE ("Dielectric coefficient of the implicit solvent");
1854 RTYPE ("gb-epsilon-solvent", ir->gb_epsilon_solvent, 80.0);
1855 CTYPE ("Salt concentration in M for Generalized Born models");
1856 RTYPE ("gb-saltconc", ir->gb_saltconc, 0.0);
1857 CTYPE ("Scaling factors used in the OBC GB model. Default values are OBC(II)");
1858 RTYPE ("gb-obc-alpha", ir->gb_obc_alpha, 1.0);
1859 RTYPE ("gb-obc-beta", ir->gb_obc_beta, 0.8);
1860 RTYPE ("gb-obc-gamma", ir->gb_obc_gamma, 4.85);
1861 RTYPE ("gb-dielectric-offset", ir->gb_dielectric_offset, 0.009);
1862 EETYPE("sa-algorithm", ir->sa_algorithm, esa_names);
1863 CTYPE ("Surface tension (kJ/mol/nm^2) for the SA (nonpolar surface) part of GBSA");
1864 CTYPE ("The value -1 will set default value for Still/HCT/OBC GB-models.");
1865 RTYPE ("sa-surface-tension", ir->sa_surface_tension, -1);
1867 /* Coupling stuff */
1868 CCTYPE ("OPTIONS FOR WEAK COUPLING ALGORITHMS");
1869 CTYPE ("Temperature coupling");
1870 EETYPE("tcoupl", ir->etc, etcoupl_names);
1871 ITYPE ("nsttcouple", ir->nsttcouple, -1);
1872 ITYPE("nh-chain-length", ir->opts.nhchainlength, NHCHAINLENGTH);
1873 EETYPE("print-nose-hoover-chain-variables", ir->bPrintNHChains, yesno_names);
1874 CTYPE ("Groups to couple separately");
1875 STYPE ("tc-grps", tcgrps, NULL);
1876 CTYPE ("Time constant (ps) and reference temperature (K)");
1877 STYPE ("tau-t", tau_t, NULL);
1878 STYPE ("ref-t", ref_t, NULL);
1879 CTYPE ("pressure coupling");
1880 EETYPE("pcoupl", ir->epc, epcoupl_names);
1881 EETYPE("pcoupltype", ir->epct, epcoupltype_names);
1882 ITYPE ("nstpcouple", ir->nstpcouple, -1);
1883 CTYPE ("Time constant (ps), compressibility (1/bar) and reference P (bar)");
1884 RTYPE ("tau-p", ir->tau_p, 1.0);
1885 STYPE ("compressibility", dumstr[0], NULL);
1886 STYPE ("ref-p", dumstr[1], NULL);
1887 CTYPE ("Scaling of reference coordinates, No, All or COM");
1888 EETYPE ("refcoord-scaling", ir->refcoord_scaling, erefscaling_names);
1891 CCTYPE ("OPTIONS FOR QMMM calculations");
1892 EETYPE("QMMM", ir->bQMMM, yesno_names);
1893 CTYPE ("Groups treated Quantum Mechanically");
1894 STYPE ("QMMM-grps", QMMM, NULL);
1895 CTYPE ("QM method");
1896 STYPE("QMmethod", QMmethod, NULL);
1897 CTYPE ("QMMM scheme");
1898 EETYPE("QMMMscheme", ir->QMMMscheme, eQMMMscheme_names);
1899 CTYPE ("QM basisset");
1900 STYPE("QMbasis", QMbasis, NULL);
1901 CTYPE ("QM charge");
1902 STYPE ("QMcharge", QMcharge, NULL);
1903 CTYPE ("QM multiplicity");
1904 STYPE ("QMmult", QMmult, NULL);
1905 CTYPE ("Surface Hopping");
1906 STYPE ("SH", bSH, NULL);
1907 CTYPE ("CAS space options");
1908 STYPE ("CASorbitals", CASorbitals, NULL);
1909 STYPE ("CASelectrons", CASelectrons, NULL);
1910 STYPE ("SAon", SAon, NULL);
1911 STYPE ("SAoff", SAoff, NULL);
1912 STYPE ("SAsteps", SAsteps, NULL);
1913 CTYPE ("Scale factor for MM charges");
1914 RTYPE ("MMChargeScaleFactor", ir->scalefactor, 1.0);
1915 CTYPE ("Optimization of QM subsystem");
1916 STYPE ("bOPT", bOPT, NULL);
1917 STYPE ("bTS", bTS, NULL);
1919 /* Simulated annealing */
1920 CCTYPE("SIMULATED ANNEALING");
1921 CTYPE ("Type of annealing for each temperature group (no/single/periodic)");
1922 STYPE ("annealing", anneal, NULL);
1923 CTYPE ("Number of time points to use for specifying annealing in each group");
1924 STYPE ("annealing-npoints", anneal_npoints, NULL);
1925 CTYPE ("List of times at the annealing points for each group");
1926 STYPE ("annealing-time", anneal_time, NULL);
1927 CTYPE ("Temp. at each annealing point, for each group.");
1928 STYPE ("annealing-temp", anneal_temp, NULL);
1931 CCTYPE ("GENERATE VELOCITIES FOR STARTUP RUN");
1932 EETYPE("gen-vel", opts->bGenVel, yesno_names);
1933 RTYPE ("gen-temp", opts->tempi, 300.0);
1934 ITYPE ("gen-seed", opts->seed, 173529);
1937 CCTYPE ("OPTIONS FOR BONDS");
1938 EETYPE("constraints", opts->nshake, constraints);
1939 CTYPE ("Type of constraint algorithm");
1940 EETYPE("constraint-algorithm", ir->eConstrAlg, econstr_names);
1941 CTYPE ("Do not constrain the start configuration");
1942 EETYPE("continuation", ir->bContinuation, yesno_names);
1943 CTYPE ("Use successive overrelaxation to reduce the number of shake iterations");
1944 EETYPE("Shake-SOR", ir->bShakeSOR, yesno_names);
1945 CTYPE ("Relative tolerance of shake");
1946 RTYPE ("shake-tol", ir->shake_tol, 0.0001);
1947 CTYPE ("Highest order in the expansion of the constraint coupling matrix");
1948 ITYPE ("lincs-order", ir->nProjOrder, 4);
1949 CTYPE ("Number of iterations in the final step of LINCS. 1 is fine for");
1950 CTYPE ("normal simulations, but use 2 to conserve energy in NVE runs.");
1951 CTYPE ("For energy minimization with constraints it should be 4 to 8.");
1952 ITYPE ("lincs-iter", ir->nLincsIter, 1);
1953 CTYPE ("Lincs will write a warning to the stderr if in one step a bond");
1954 CTYPE ("rotates over more degrees than");
1955 RTYPE ("lincs-warnangle", ir->LincsWarnAngle, 30.0);
1956 CTYPE ("Convert harmonic bonds to morse potentials");
1957 EETYPE("morse", opts->bMorse, yesno_names);
1959 /* Energy group exclusions */
1960 CCTYPE ("ENERGY GROUP EXCLUSIONS");
1961 CTYPE ("Pairs of energy groups for which all non-bonded interactions are excluded");
1962 STYPE ("energygrp-excl", egpexcl, NULL);
1966 CTYPE ("Number of walls, type, atom types, densities and box-z scale factor for Ewald");
1967 ITYPE ("nwall", ir->nwall, 0);
1968 EETYPE("wall-type", ir->wall_type, ewt_names);
1969 RTYPE ("wall-r-linpot", ir->wall_r_linpot, -1);
1970 STYPE ("wall-atomtype", wall_atomtype, NULL);
1971 STYPE ("wall-density", wall_density, NULL);
1972 RTYPE ("wall-ewald-zfac", ir->wall_ewald_zfac, 3);
1975 CCTYPE("COM PULLING");
1976 CTYPE("Pull type: no, umbrella, constraint or constant-force");
1977 EETYPE("pull", ir->ePull, epull_names);
1978 if (ir->ePull != epullNO)
1981 pull_grp = read_pullparams(&ninp, &inp, ir->pull, &opts->pull_start, wi);
1984 /* Enforced rotation */
1985 CCTYPE("ENFORCED ROTATION");
1986 CTYPE("Enforced rotation: No or Yes");
1987 EETYPE("rotation", ir->bRot, yesno_names);
1991 rot_grp = read_rotparams(&ninp, &inp, ir->rot, wi);
1995 CCTYPE("NMR refinement stuff");
1996 CTYPE ("Distance restraints type: No, Simple or Ensemble");
1997 EETYPE("disre", ir->eDisre, edisre_names);
1998 CTYPE ("Force weighting of pairs in one distance restraint: Conservative or Equal");
1999 EETYPE("disre-weighting", ir->eDisreWeighting, edisreweighting_names);
2000 CTYPE ("Use sqrt of the time averaged times the instantaneous violation");
2001 EETYPE("disre-mixed", ir->bDisreMixed, yesno_names);
2002 RTYPE ("disre-fc", ir->dr_fc, 1000.0);
2003 RTYPE ("disre-tau", ir->dr_tau, 0.0);
2004 CTYPE ("Output frequency for pair distances to energy file");
2005 ITYPE ("nstdisreout", ir->nstdisreout, 100);
2006 CTYPE ("Orientation restraints: No or Yes");
2007 EETYPE("orire", opts->bOrire, yesno_names);
2008 CTYPE ("Orientation restraints force constant and tau for time averaging");
2009 RTYPE ("orire-fc", ir->orires_fc, 0.0);
2010 RTYPE ("orire-tau", ir->orires_tau, 0.0);
2011 STYPE ("orire-fitgrp", orirefitgrp, NULL);
2012 CTYPE ("Output frequency for trace(SD) and S to energy file");
2013 ITYPE ("nstorireout", ir->nstorireout, 100);
2015 /* free energy variables */
2016 CCTYPE ("Free energy variables");
2017 EETYPE("free-energy", ir->efep, efep_names);
2018 STYPE ("couple-moltype", couple_moltype, NULL);
2019 EETYPE("couple-lambda0", opts->couple_lam0, couple_lam);
2020 EETYPE("couple-lambda1", opts->couple_lam1, couple_lam);
2021 EETYPE("couple-intramol", opts->bCoupleIntra, yesno_names);
2023 RTYPE ("init-lambda", fep->init_lambda, -1); /* start with -1 so
2025 it was not entered */
2026 ITYPE ("init-lambda-state", fep->init_fep_state, -1);
2027 RTYPE ("delta-lambda", fep->delta_lambda, 0.0);
2028 ITYPE ("nstdhdl", fep->nstdhdl, 50);
2029 STYPE ("fep-lambdas", fep_lambda[efptFEP], NULL);
2030 STYPE ("mass-lambdas", fep_lambda[efptMASS], NULL);
2031 STYPE ("coul-lambdas", fep_lambda[efptCOUL], NULL);
2032 STYPE ("vdw-lambdas", fep_lambda[efptVDW], NULL);
2033 STYPE ("bonded-lambdas", fep_lambda[efptBONDED], NULL);
2034 STYPE ("restraint-lambdas", fep_lambda[efptRESTRAINT], NULL);
2035 STYPE ("temperature-lambdas", fep_lambda[efptTEMPERATURE], NULL);
2036 ITYPE ("calc-lambda-neighbors", fep->lambda_neighbors, 1);
2037 STYPE ("init-lambda-weights", lambda_weights, NULL);
2038 EETYPE("dhdl-print-energy", fep->bPrintEnergy, yesno_names);
2039 RTYPE ("sc-alpha", fep->sc_alpha, 0.0);
2040 ITYPE ("sc-power", fep->sc_power, 1);
2041 RTYPE ("sc-r-power", fep->sc_r_power, 6.0);
2042 RTYPE ("sc-sigma", fep->sc_sigma, 0.3);
2043 EETYPE("sc-coul", fep->bScCoul, yesno_names);
2044 ITYPE ("dh_hist_size", fep->dh_hist_size, 0);
2045 RTYPE ("dh_hist_spacing", fep->dh_hist_spacing, 0.1);
2046 EETYPE("separate-dhdl-file", fep->separate_dhdl_file,
2047 separate_dhdl_file_names);
2048 EETYPE("dhdl-derivatives", fep->dhdl_derivatives, dhdl_derivatives_names);
2049 ITYPE ("dh_hist_size", fep->dh_hist_size, 0);
2050 RTYPE ("dh_hist_spacing", fep->dh_hist_spacing, 0.1);
2052 /* Non-equilibrium MD stuff */
2053 CCTYPE("Non-equilibrium MD stuff");
2054 STYPE ("acc-grps", accgrps, NULL);
2055 STYPE ("accelerate", acc, NULL);
2056 STYPE ("freezegrps", freeze, NULL);
2057 STYPE ("freezedim", frdim, NULL);
2058 RTYPE ("cos-acceleration", ir->cos_accel, 0);
2059 STYPE ("deform", deform, NULL);
2061 /* simulated tempering variables */
2062 CCTYPE("simulated tempering variables");
2063 EETYPE("simulated-tempering", ir->bSimTemp, yesno_names);
2064 EETYPE("simulated-tempering-scaling", ir->simtempvals->eSimTempScale, esimtemp_names);
2065 RTYPE("sim-temp-low", ir->simtempvals->simtemp_low, 300.0);
2066 RTYPE("sim-temp-high", ir->simtempvals->simtemp_high, 300.0);
2068 /* expanded ensemble variables */
2069 if (ir->efep == efepEXPANDED || ir->bSimTemp)
2071 read_expandedparams(&ninp, &inp, expand, wi);
2074 /* Electric fields */
2075 CCTYPE("Electric fields");
2076 CTYPE ("Format is number of terms (int) and for all terms an amplitude (real)");
2077 CTYPE ("and a phase angle (real)");
2078 STYPE ("E-x", efield_x, NULL);
2079 STYPE ("E-xt", efield_xt, NULL);
2080 STYPE ("E-y", efield_y, NULL);
2081 STYPE ("E-yt", efield_yt, NULL);
2082 STYPE ("E-z", efield_z, NULL);
2083 STYPE ("E-zt", efield_zt, NULL);
2085 /* AdResS defined thingies */
2086 CCTYPE ("AdResS parameters");
2087 EETYPE("adress", ir->bAdress, yesno_names);
2090 snew(ir->adress, 1);
2091 read_adressparams(&ninp, &inp, ir->adress, wi);
2094 /* User defined thingies */
2095 CCTYPE ("User defined thingies");
2096 STYPE ("user1-grps", user1, NULL);
2097 STYPE ("user2-grps", user2, NULL);
2098 ITYPE ("userint1", ir->userint1, 0);
2099 ITYPE ("userint2", ir->userint2, 0);
2100 ITYPE ("userint3", ir->userint3, 0);
2101 ITYPE ("userint4", ir->userint4, 0);
2102 RTYPE ("userreal1", ir->userreal1, 0);
2103 RTYPE ("userreal2", ir->userreal2, 0);
2104 RTYPE ("userreal3", ir->userreal3, 0);
2105 RTYPE ("userreal4", ir->userreal4, 0);
2108 write_inpfile(mdparout, ninp, inp, FALSE, wi);
2109 for (i = 0; (i < ninp); i++)
2112 sfree(inp[i].value);
2116 /* Process options if necessary */
2117 for (m = 0; m < 2; m++)
2119 for (i = 0; i < 2*DIM; i++)
2128 if (sscanf(dumstr[m], "%lf", &(dumdub[m][XX])) != 1)
2130 warning_error(wi, "Pressure coupling not enough values (I need 1)");
2132 dumdub[m][YY] = dumdub[m][ZZ] = dumdub[m][XX];
2134 case epctSEMIISOTROPIC:
2135 case epctSURFACETENSION:
2136 if (sscanf(dumstr[m], "%lf%lf",
2137 &(dumdub[m][XX]), &(dumdub[m][ZZ])) != 2)
2139 warning_error(wi, "Pressure coupling not enough values (I need 2)");
2141 dumdub[m][YY] = dumdub[m][XX];
2143 case epctANISOTROPIC:
2144 if (sscanf(dumstr[m], "%lf%lf%lf%lf%lf%lf",
2145 &(dumdub[m][XX]), &(dumdub[m][YY]), &(dumdub[m][ZZ]),
2146 &(dumdub[m][3]), &(dumdub[m][4]), &(dumdub[m][5])) != 6)
2148 warning_error(wi, "Pressure coupling not enough values (I need 6)");
2152 gmx_fatal(FARGS, "Pressure coupling type %s not implemented yet",
2153 epcoupltype_names[ir->epct]);
2157 clear_mat(ir->ref_p);
2158 clear_mat(ir->compress);
2159 for (i = 0; i < DIM; i++)
2161 ir->ref_p[i][i] = dumdub[1][i];
2162 ir->compress[i][i] = dumdub[0][i];
2164 if (ir->epct == epctANISOTROPIC)
2166 ir->ref_p[XX][YY] = dumdub[1][3];
2167 ir->ref_p[XX][ZZ] = dumdub[1][4];
2168 ir->ref_p[YY][ZZ] = dumdub[1][5];
2169 if (ir->ref_p[XX][YY] != 0 && ir->ref_p[XX][ZZ] != 0 && ir->ref_p[YY][ZZ] != 0)
2171 warning(wi, "All off-diagonal reference pressures are non-zero. Are you sure you want to apply a threefold shear stress?\n");
2173 ir->compress[XX][YY] = dumdub[0][3];
2174 ir->compress[XX][ZZ] = dumdub[0][4];
2175 ir->compress[YY][ZZ] = dumdub[0][5];
2176 for (i = 0; i < DIM; i++)
2178 for (m = 0; m < i; m++)
2180 ir->ref_p[i][m] = ir->ref_p[m][i];
2181 ir->compress[i][m] = ir->compress[m][i];
2186 if (ir->comm_mode == ecmNO)
2191 opts->couple_moltype = NULL;
2192 if (strlen(couple_moltype) > 0)
2194 if (ir->efep != efepNO)
2196 opts->couple_moltype = strdup(couple_moltype);
2197 if (opts->couple_lam0 == opts->couple_lam1)
2199 warning(wi, "The lambda=0 and lambda=1 states for coupling are identical");
2201 if (ir->eI == eiMD && (opts->couple_lam0 == ecouplamNONE ||
2202 opts->couple_lam1 == ecouplamNONE))
2204 warning(wi, "For proper sampling of the (nearly) decoupled state, stochastic dynamics should be used");
2209 warning(wi, "Can not couple a molecule with free_energy = no");
2212 /* FREE ENERGY AND EXPANDED ENSEMBLE OPTIONS */
2213 if (ir->efep != efepNO)
2215 if (fep->delta_lambda > 0)
2217 ir->efep = efepSLOWGROWTH;
2223 fep->bPrintEnergy = TRUE;
2224 /* always print out the energy to dhdl if we are doing expanded ensemble, since we need the total energy
2225 if the temperature is changing. */
2228 if ((ir->efep != efepNO) || ir->bSimTemp)
2230 ir->bExpanded = FALSE;
2231 if ((ir->efep == efepEXPANDED) || ir->bSimTemp)
2233 ir->bExpanded = TRUE;
2235 do_fep_params(ir, fep_lambda, lambda_weights);
2236 if (ir->bSimTemp) /* done after fep params */
2238 do_simtemp_params(ir);
2243 ir->fepvals->n_lambda = 0;
2246 /* WALL PARAMETERS */
2248 do_wall_params(ir, wall_atomtype, wall_density, opts);
2250 /* ORIENTATION RESTRAINT PARAMETERS */
2252 if (opts->bOrire && str_nelem(orirefitgrp, MAXPTR, NULL) != 1)
2254 warning_error(wi, "ERROR: Need one orientation restraint fit group\n");
2257 /* DEFORMATION PARAMETERS */
2259 clear_mat(ir->deform);
2260 for (i = 0; i < 6; i++)
2264 m = sscanf(deform, "%lf %lf %lf %lf %lf %lf",
2265 &(dumdub[0][0]), &(dumdub[0][1]), &(dumdub[0][2]),
2266 &(dumdub[0][3]), &(dumdub[0][4]), &(dumdub[0][5]));
2267 for (i = 0; i < 3; i++)
2269 ir->deform[i][i] = dumdub[0][i];
2271 ir->deform[YY][XX] = dumdub[0][3];
2272 ir->deform[ZZ][XX] = dumdub[0][4];
2273 ir->deform[ZZ][YY] = dumdub[0][5];
2274 if (ir->epc != epcNO)
2276 for (i = 0; i < 3; i++)
2278 for (j = 0; j <= i; j++)
2280 if (ir->deform[i][j] != 0 && ir->compress[i][j] != 0)
2282 warning_error(wi, "A box element has deform set and compressibility > 0");
2286 for (i = 0; i < 3; i++)
2288 for (j = 0; j < i; j++)
2290 if (ir->deform[i][j] != 0)
2292 for (m = j; m < DIM; m++)
2294 if (ir->compress[m][j] != 0)
2296 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.");
2297 warning(wi, warn_buf);
2309 static int search_QMstring(char *s, int ng, const char *gn[])
2311 /* same as normal search_string, but this one searches QM strings */
2314 for (i = 0; (i < ng); i++)
2316 if (gmx_strcasecmp(s, gn[i]) == 0)
2322 gmx_fatal(FARGS, "this QM method or basisset (%s) is not implemented\n!", s);
2326 } /* search_QMstring */
2329 int search_string(char *s, int ng, char *gn[])
2333 for (i = 0; (i < ng); i++)
2335 if (gmx_strcasecmp(s, gn[i]) == 0)
2342 "Group %s referenced in the .mdp file was not found in the index file.\n"
2343 "Group names must match either [moleculetype] names or custom index group\n"
2344 "names, in which case you must supply an index file to the '-n' option\n"
2351 static gmx_bool do_numbering(int natoms, gmx_groups_t *groups, int ng, char *ptrs[],
2352 t_blocka *block, char *gnames[],
2353 int gtype, int restnm,
2354 int grptp, gmx_bool bVerbose,
2357 unsigned short *cbuf;
2358 t_grps *grps = &(groups->grps[gtype]);
2359 int i, j, gid, aj, ognr, ntot = 0;
2362 char warn_buf[STRLEN];
2366 fprintf(debug, "Starting numbering %d groups of type %d\n", ng, gtype);
2369 title = gtypes[gtype];
2372 /* Mark all id's as not set */
2373 for (i = 0; (i < natoms); i++)
2378 snew(grps->nm_ind, ng+1); /* +1 for possible rest group */
2379 for (i = 0; (i < ng); i++)
2381 /* Lookup the group name in the block structure */
2382 gid = search_string(ptrs[i], block->nr, gnames);
2383 if ((grptp != egrptpONE) || (i == 0))
2385 grps->nm_ind[grps->nr++] = gid;
2389 fprintf(debug, "Found gid %d for group %s\n", gid, ptrs[i]);
2392 /* Now go over the atoms in the group */
2393 for (j = block->index[gid]; (j < block->index[gid+1]); j++)
2398 /* Range checking */
2399 if ((aj < 0) || (aj >= natoms))
2401 gmx_fatal(FARGS, "Invalid atom number %d in indexfile", aj);
2403 /* Lookup up the old group number */
2407 gmx_fatal(FARGS, "Atom %d in multiple %s groups (%d and %d)",
2408 aj+1, title, ognr+1, i+1);
2412 /* Store the group number in buffer */
2413 if (grptp == egrptpONE)
2426 /* Now check whether we have done all atoms */
2430 if (grptp == egrptpALL)
2432 gmx_fatal(FARGS, "%d atoms are not part of any of the %s groups",
2433 natoms-ntot, title);
2435 else if (grptp == egrptpPART)
2437 sprintf(warn_buf, "%d atoms are not part of any of the %s groups",
2438 natoms-ntot, title);
2439 warning_note(wi, warn_buf);
2441 /* Assign all atoms currently unassigned to a rest group */
2442 for (j = 0; (j < natoms); j++)
2444 if (cbuf[j] == NOGID)
2450 if (grptp != egrptpPART)
2455 "Making dummy/rest group for %s containing %d elements\n",
2456 title, natoms-ntot);
2458 /* Add group name "rest" */
2459 grps->nm_ind[grps->nr] = restnm;
2461 /* Assign the rest name to all atoms not currently assigned to a group */
2462 for (j = 0; (j < natoms); j++)
2464 if (cbuf[j] == NOGID)
2473 if (grps->nr == 1 && (ntot == 0 || ntot == natoms))
2475 /* All atoms are part of one (or no) group, no index required */
2476 groups->ngrpnr[gtype] = 0;
2477 groups->grpnr[gtype] = NULL;
2481 groups->ngrpnr[gtype] = natoms;
2482 snew(groups->grpnr[gtype], natoms);
2483 for (j = 0; (j < natoms); j++)
2485 groups->grpnr[gtype][j] = cbuf[j];
2491 return (bRest && grptp == egrptpPART);
2494 static void calc_nrdf(gmx_mtop_t *mtop, t_inputrec *ir, char **gnames)
2497 gmx_groups_t *groups;
2499 int natoms, ai, aj, i, j, d, g, imin, jmin, nc;
2501 int *nrdf2, *na_vcm, na_tot;
2502 double *nrdf_tc, *nrdf_vcm, nrdf_uc, n_sub = 0;
2503 gmx_mtop_atomloop_all_t aloop;
2505 int mb, mol, ftype, as;
2506 gmx_molblock_t *molb;
2507 gmx_moltype_t *molt;
2510 * First calc 3xnr-atoms for each group
2511 * then subtract half a degree of freedom for each constraint
2513 * Only atoms and nuclei contribute to the degrees of freedom...
2518 groups = &mtop->groups;
2519 natoms = mtop->natoms;
2521 /* Allocate one more for a possible rest group */
2522 /* We need to sum degrees of freedom into doubles,
2523 * since floats give too low nrdf's above 3 million atoms.
2525 snew(nrdf_tc, groups->grps[egcTC].nr+1);
2526 snew(nrdf_vcm, groups->grps[egcVCM].nr+1);
2527 snew(na_vcm, groups->grps[egcVCM].nr+1);
2529 for (i = 0; i < groups->grps[egcTC].nr; i++)
2533 for (i = 0; i < groups->grps[egcVCM].nr+1; i++)
2538 snew(nrdf2, natoms);
2539 aloop = gmx_mtop_atomloop_all_init(mtop);
2540 while (gmx_mtop_atomloop_all_next(aloop, &i, &atom))
2543 if (atom->ptype == eptAtom || atom->ptype == eptNucleus)
2545 g = ggrpnr(groups, egcFREEZE, i);
2546 /* Double count nrdf for particle i */
2547 for (d = 0; d < DIM; d++)
2549 if (opts->nFreeze[g][d] == 0)
2554 nrdf_tc [ggrpnr(groups, egcTC, i)] += 0.5*nrdf2[i];
2555 nrdf_vcm[ggrpnr(groups, egcVCM, i)] += 0.5*nrdf2[i];
2560 for (mb = 0; mb < mtop->nmolblock; mb++)
2562 molb = &mtop->molblock[mb];
2563 molt = &mtop->moltype[molb->type];
2564 atom = molt->atoms.atom;
2565 for (mol = 0; mol < molb->nmol; mol++)
2567 for (ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++)
2569 ia = molt->ilist[ftype].iatoms;
2570 for (i = 0; i < molt->ilist[ftype].nr; )
2572 /* Subtract degrees of freedom for the constraints,
2573 * if the particles still have degrees of freedom left.
2574 * If one of the particles is a vsite or a shell, then all
2575 * constraint motion will go there, but since they do not
2576 * contribute to the constraints the degrees of freedom do not
2581 if (((atom[ia[1]].ptype == eptNucleus) ||
2582 (atom[ia[1]].ptype == eptAtom)) &&
2583 ((atom[ia[2]].ptype == eptNucleus) ||
2584 (atom[ia[2]].ptype == eptAtom)))
2602 imin = min(imin, nrdf2[ai]);
2603 jmin = min(jmin, nrdf2[aj]);
2606 nrdf_tc [ggrpnr(groups, egcTC, ai)] -= 0.5*imin;
2607 nrdf_tc [ggrpnr(groups, egcTC, aj)] -= 0.5*jmin;
2608 nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5*imin;
2609 nrdf_vcm[ggrpnr(groups, egcVCM, aj)] -= 0.5*jmin;
2611 ia += interaction_function[ftype].nratoms+1;
2612 i += interaction_function[ftype].nratoms+1;
2615 ia = molt->ilist[F_SETTLE].iatoms;
2616 for (i = 0; i < molt->ilist[F_SETTLE].nr; )
2618 /* Subtract 1 dof from every atom in the SETTLE */
2619 for (j = 0; j < 3; j++)
2622 imin = min(2, nrdf2[ai]);
2624 nrdf_tc [ggrpnr(groups, egcTC, ai)] -= 0.5*imin;
2625 nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5*imin;
2630 as += molt->atoms.nr;
2634 if (ir->ePull == epullCONSTRAINT)
2636 /* Correct nrdf for the COM constraints.
2637 * We correct using the TC and VCM group of the first atom
2638 * in the reference and pull group. If atoms in one pull group
2639 * belong to different TC or VCM groups it is anyhow difficult
2640 * to determine the optimal nrdf assignment.
2643 if (pull->eGeom == epullgPOS)
2646 for (i = 0; i < DIM; i++)
2658 for (i = 0; i < pull->ngrp; i++)
2661 if (pull->grp[0].nat > 0)
2663 /* Subtract 1/2 dof from the reference group */
2664 ai = pull->grp[0].ind[0];
2665 if (nrdf_tc[ggrpnr(groups, egcTC, ai)] > 1)
2667 nrdf_tc [ggrpnr(groups, egcTC, ai)] -= 0.5;
2668 nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5;
2672 /* Subtract 1/2 dof from the pulled group */
2673 ai = pull->grp[1+i].ind[0];
2674 nrdf_tc [ggrpnr(groups, egcTC, ai)] -= 0.5*imin;
2675 nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5*imin;
2676 if (nrdf_tc[ggrpnr(groups, egcTC, ai)] < 0)
2678 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)]]);
2683 if (ir->nstcomm != 0)
2685 /* Subtract 3 from the number of degrees of freedom in each vcm group
2686 * when com translation is removed and 6 when rotation is removed
2689 switch (ir->comm_mode)
2692 n_sub = ndof_com(ir);
2699 gmx_incons("Checking comm_mode");
2702 for (i = 0; i < groups->grps[egcTC].nr; i++)
2704 /* Count the number of atoms of TC group i for every VCM group */
2705 for (j = 0; j < groups->grps[egcVCM].nr+1; j++)
2710 for (ai = 0; ai < natoms; ai++)
2712 if (ggrpnr(groups, egcTC, ai) == i)
2714 na_vcm[ggrpnr(groups, egcVCM, ai)]++;
2718 /* Correct for VCM removal according to the fraction of each VCM
2719 * group present in this TC group.
2721 nrdf_uc = nrdf_tc[i];
2724 fprintf(debug, "T-group[%d] nrdf_uc = %g, n_sub = %g\n",
2728 for (j = 0; j < groups->grps[egcVCM].nr+1; j++)
2730 if (nrdf_vcm[j] > n_sub)
2732 nrdf_tc[i] += nrdf_uc*((double)na_vcm[j]/(double)na_tot)*
2733 (nrdf_vcm[j] - n_sub)/nrdf_vcm[j];
2737 fprintf(debug, " nrdf_vcm[%d] = %g, nrdf = %g\n",
2738 j, nrdf_vcm[j], nrdf_tc[i]);
2743 for (i = 0; (i < groups->grps[egcTC].nr); i++)
2745 opts->nrdf[i] = nrdf_tc[i];
2746 if (opts->nrdf[i] < 0)
2751 "Number of degrees of freedom in T-Coupling group %s is %.2f\n",
2752 gnames[groups->grps[egcTC].nm_ind[i]], opts->nrdf[i]);
2761 static void decode_cos(char *s, t_cosines *cosine)
2764 char format[STRLEN], f1[STRLEN];
2776 sscanf(t, "%d", &(cosine->n));
2783 snew(cosine->a, cosine->n);
2784 snew(cosine->phi, cosine->n);
2786 sprintf(format, "%%*d");
2787 for (i = 0; (i < cosine->n); i++)
2790 strcat(f1, "%lf%lf");
2791 if (sscanf(t, f1, &a, &phi) < 2)
2793 gmx_fatal(FARGS, "Invalid input for electric field shift: '%s'", t);
2796 cosine->phi[i] = phi;
2797 strcat(format, "%*lf%*lf");
2804 static gmx_bool do_egp_flag(t_inputrec *ir, gmx_groups_t *groups,
2805 const char *option, const char *val, int flag)
2807 /* The maximum number of energy group pairs would be MAXPTR*(MAXPTR+1)/2.
2808 * But since this is much larger than STRLEN, such a line can not be parsed.
2809 * The real maximum is the number of names that fit in a string: STRLEN/2.
2811 #define EGP_MAX (STRLEN/2)
2812 int nelem, i, j, k, nr;
2813 char *names[EGP_MAX];
2817 gnames = groups->grpname;
2819 nelem = str_nelem(val, EGP_MAX, names);
2822 gmx_fatal(FARGS, "The number of groups for %s is odd", option);
2824 nr = groups->grps[egcENER].nr;
2826 for (i = 0; i < nelem/2; i++)
2830 gmx_strcasecmp(names[2*i], *(gnames[groups->grps[egcENER].nm_ind[j]])))
2836 gmx_fatal(FARGS, "%s in %s is not an energy group\n",
2837 names[2*i], option);
2841 gmx_strcasecmp(names[2*i+1], *(gnames[groups->grps[egcENER].nm_ind[k]])))
2847 gmx_fatal(FARGS, "%s in %s is not an energy group\n",
2848 names[2*i+1], option);
2850 if ((j < nr) && (k < nr))
2852 ir->opts.egp_flags[nr*j+k] |= flag;
2853 ir->opts.egp_flags[nr*k+j] |= flag;
2861 void do_index(const char* mdparin, const char *ndx,
2864 t_inputrec *ir, rvec *v,
2868 gmx_groups_t *groups;
2872 char warnbuf[STRLEN], **gnames;
2873 int nr, ntcg, ntau_t, nref_t, nacc, nofg, nSA, nSA_points, nSA_time, nSA_temp;
2876 int nacg, nfreeze, nfrdim, nenergy, nvcm, nuser;
2877 char *ptr1[MAXPTR], *ptr2[MAXPTR], *ptr3[MAXPTR];
2878 int i, j, k, restnm;
2880 gmx_bool bExcl, bTable, bSetTCpar, bAnneal, bRest;
2881 int nQMmethod, nQMbasis, nQMcharge, nQMmult, nbSH, nCASorb, nCASelec,
2882 nSAon, nSAoff, nSAsteps, nQMg, nbOPT, nbTS;
2883 char warn_buf[STRLEN];
2887 fprintf(stderr, "processing index file...\n");
2893 snew(grps->index, 1);
2895 atoms_all = gmx_mtop_global_atoms(mtop);
2896 analyse(&atoms_all, grps, &gnames, FALSE, TRUE);
2897 free_t_atoms(&atoms_all, FALSE);
2901 grps = init_index(ndx, &gnames);
2904 groups = &mtop->groups;
2905 natoms = mtop->natoms;
2906 symtab = &mtop->symtab;
2908 snew(groups->grpname, grps->nr+1);
2910 for (i = 0; (i < grps->nr); i++)
2912 groups->grpname[i] = put_symtab(symtab, gnames[i]);
2914 groups->grpname[i] = put_symtab(symtab, "rest");
2916 srenew(gnames, grps->nr+1);
2917 gnames[restnm] = *(groups->grpname[i]);
2918 groups->ngrpname = grps->nr+1;
2920 set_warning_line(wi, mdparin, -1);
2922 ntau_t = str_nelem(tau_t, MAXPTR, ptr1);
2923 nref_t = str_nelem(ref_t, MAXPTR, ptr2);
2924 ntcg = str_nelem(tcgrps, MAXPTR, ptr3);
2925 if ((ntau_t != ntcg) || (nref_t != ntcg))
2927 gmx_fatal(FARGS, "Invalid T coupling input: %d groups, %d ref-t values and "
2928 "%d tau-t values", ntcg, nref_t, ntau_t);
2931 bSetTCpar = (ir->etc || EI_SD(ir->eI) || ir->eI == eiBD || EI_TPI(ir->eI));
2932 do_numbering(natoms, groups, ntcg, ptr3, grps, gnames, egcTC,
2933 restnm, bSetTCpar ? egrptpALL : egrptpALL_GENREST, bVerbose, wi);
2934 nr = groups->grps[egcTC].nr;
2936 snew(ir->opts.nrdf, nr);
2937 snew(ir->opts.tau_t, nr);
2938 snew(ir->opts.ref_t, nr);
2939 if (ir->eI == eiBD && ir->bd_fric == 0)
2941 fprintf(stderr, "bd-fric=0, so tau-t will be used as the inverse friction constant(s)\n");
2948 gmx_fatal(FARGS, "Not enough ref-t and tau-t values!");
2952 for (i = 0; (i < nr); i++)
2954 ir->opts.tau_t[i] = strtod(ptr1[i], NULL);
2955 if ((ir->eI == eiBD || ir->eI == eiSD2) && ir->opts.tau_t[i] <= 0)
2957 sprintf(warn_buf, "With integrator %s tau-t should be larger than 0", ei_names[ir->eI]);
2958 warning_error(wi, warn_buf);
2961 if (ir->etc != etcVRESCALE && ir->opts.tau_t[i] == 0)
2963 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.");
2966 if (ir->opts.tau_t[i] >= 0)
2968 tau_min = min(tau_min, ir->opts.tau_t[i]);
2971 if (ir->etc != etcNO && ir->nsttcouple == -1)
2973 ir->nsttcouple = ir_optimal_nsttcouple(ir);
2978 if ((ir->etc == etcNOSEHOOVER) && (ir->epc == epcBERENDSEN))
2980 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");
2982 if ((ir->epc == epcMTTK) && (ir->etc > etcNO))
2984 if (ir->nstpcouple != ir->nsttcouple)
2986 int mincouple = min(ir->nstpcouple, ir->nsttcouple);
2987 ir->nstpcouple = ir->nsttcouple = mincouple;
2988 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);
2989 warning_note(wi, warn_buf);
2993 /* velocity verlet with averaged kinetic energy KE = 0.5*(v(t+1/2) - v(t-1/2)) is implemented
2994 primarily for testing purposes, and does not work with temperature coupling other than 1 */
2996 if (ETC_ANDERSEN(ir->etc))
2998 if (ir->nsttcouple != 1)
3001 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");
3002 warning_note(wi, warn_buf);
3005 nstcmin = tcouple_min_integration_steps(ir->etc);
3008 if (tau_min/(ir->delta_t*ir->nsttcouple) < nstcmin)
3010 sprintf(warn_buf, "For proper integration of the %s thermostat, tau-t (%g) should be at least %d times larger than nsttcouple*dt (%g)",
3011 ETCOUPLTYPE(ir->etc),
3013 ir->nsttcouple*ir->delta_t);
3014 warning(wi, warn_buf);
3017 for (i = 0; (i < nr); i++)
3019 ir->opts.ref_t[i] = strtod(ptr2[i], NULL);
3020 if (ir->opts.ref_t[i] < 0)
3022 gmx_fatal(FARGS, "ref-t for group %d negative", i);
3025 /* set the lambda mc temperature to the md integrator temperature (which should be defined
3026 if we are in this conditional) if mc_temp is negative */
3027 if (ir->expandedvals->mc_temp < 0)
3029 ir->expandedvals->mc_temp = ir->opts.ref_t[0]; /*for now, set to the first reft */
3033 /* Simulated annealing for each group. There are nr groups */
3034 nSA = str_nelem(anneal, MAXPTR, ptr1);
3035 if (nSA == 1 && (ptr1[0][0] == 'n' || ptr1[0][0] == 'N'))
3039 if (nSA > 0 && nSA != nr)
3041 gmx_fatal(FARGS, "Not enough annealing values: %d (for %d groups)\n", nSA, nr);
3045 snew(ir->opts.annealing, nr);
3046 snew(ir->opts.anneal_npoints, nr);
3047 snew(ir->opts.anneal_time, nr);
3048 snew(ir->opts.anneal_temp, nr);
3049 for (i = 0; i < nr; i++)
3051 ir->opts.annealing[i] = eannNO;
3052 ir->opts.anneal_npoints[i] = 0;
3053 ir->opts.anneal_time[i] = NULL;
3054 ir->opts.anneal_temp[i] = NULL;
3059 for (i = 0; i < nr; i++)
3061 if (ptr1[i][0] == 'n' || ptr1[i][0] == 'N')
3063 ir->opts.annealing[i] = eannNO;
3065 else if (ptr1[i][0] == 's' || ptr1[i][0] == 'S')
3067 ir->opts.annealing[i] = eannSINGLE;
3070 else if (ptr1[i][0] == 'p' || ptr1[i][0] == 'P')
3072 ir->opts.annealing[i] = eannPERIODIC;
3078 /* Read the other fields too */
3079 nSA_points = str_nelem(anneal_npoints, MAXPTR, ptr1);
3080 if (nSA_points != nSA)
3082 gmx_fatal(FARGS, "Found %d annealing-npoints values for %d groups\n", nSA_points, nSA);
3084 for (k = 0, i = 0; i < nr; i++)
3086 ir->opts.anneal_npoints[i] = strtol(ptr1[i], NULL, 10);
3087 if (ir->opts.anneal_npoints[i] == 1)
3089 gmx_fatal(FARGS, "Please specify at least a start and an end point for annealing\n");
3091 snew(ir->opts.anneal_time[i], ir->opts.anneal_npoints[i]);
3092 snew(ir->opts.anneal_temp[i], ir->opts.anneal_npoints[i]);
3093 k += ir->opts.anneal_npoints[i];
3096 nSA_time = str_nelem(anneal_time, MAXPTR, ptr1);
3099 gmx_fatal(FARGS, "Found %d annealing-time values, wanter %d\n", nSA_time, k);
3101 nSA_temp = str_nelem(anneal_temp, MAXPTR, ptr2);
3104 gmx_fatal(FARGS, "Found %d annealing-temp values, wanted %d\n", nSA_temp, k);
3107 for (i = 0, k = 0; i < nr; i++)
3110 for (j = 0; j < ir->opts.anneal_npoints[i]; j++)
3112 ir->opts.anneal_time[i][j] = strtod(ptr1[k], NULL);
3113 ir->opts.anneal_temp[i][j] = strtod(ptr2[k], NULL);
3116 if (ir->opts.anneal_time[i][0] > (ir->init_t+GMX_REAL_EPS))
3118 gmx_fatal(FARGS, "First time point for annealing > init_t.\n");
3124 if (ir->opts.anneal_time[i][j] < ir->opts.anneal_time[i][j-1])
3126 gmx_fatal(FARGS, "Annealing timepoints out of order: t=%f comes after t=%f\n",
3127 ir->opts.anneal_time[i][j], ir->opts.anneal_time[i][j-1]);
3130 if (ir->opts.anneal_temp[i][j] < 0)
3132 gmx_fatal(FARGS, "Found negative temperature in annealing: %f\n", ir->opts.anneal_temp[i][j]);
3137 /* Print out some summary information, to make sure we got it right */
3138 for (i = 0, k = 0; i < nr; i++)
3140 if (ir->opts.annealing[i] != eannNO)
3142 j = groups->grps[egcTC].nm_ind[i];
3143 fprintf(stderr, "Simulated annealing for group %s: %s, %d timepoints\n",
3144 *(groups->grpname[j]), eann_names[ir->opts.annealing[i]],
3145 ir->opts.anneal_npoints[i]);
3146 fprintf(stderr, "Time (ps) Temperature (K)\n");
3147 /* All terms except the last one */
3148 for (j = 0; j < (ir->opts.anneal_npoints[i]-1); j++)
3150 fprintf(stderr, "%9.1f %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3153 /* Finally the last one */
3154 j = ir->opts.anneal_npoints[i]-1;
3155 if (ir->opts.annealing[i] == eannSINGLE)
3157 fprintf(stderr, "%9.1f- %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3161 fprintf(stderr, "%9.1f %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3162 if (fabs(ir->opts.anneal_temp[i][j]-ir->opts.anneal_temp[i][0]) > GMX_REAL_EPS)
3164 warning_note(wi, "There is a temperature jump when your annealing loops back.\n");
3173 if (ir->ePull != epullNO)
3175 make_pull_groups(ir->pull, pull_grp, grps, gnames);
3180 make_rotation_groups(ir->rot, rot_grp, grps, gnames);
3183 nacc = str_nelem(acc, MAXPTR, ptr1);
3184 nacg = str_nelem(accgrps, MAXPTR, ptr2);
3185 if (nacg*DIM != nacc)
3187 gmx_fatal(FARGS, "Invalid Acceleration input: %d groups and %d acc. values",
3190 do_numbering(natoms, groups, nacg, ptr2, grps, gnames, egcACC,
3191 restnm, egrptpALL_GENREST, bVerbose, wi);
3192 nr = groups->grps[egcACC].nr;
3193 snew(ir->opts.acc, nr);
3194 ir->opts.ngacc = nr;
3196 for (i = k = 0; (i < nacg); i++)
3198 for (j = 0; (j < DIM); j++, k++)
3200 ir->opts.acc[i][j] = strtod(ptr1[k], NULL);
3203 for (; (i < nr); i++)
3205 for (j = 0; (j < DIM); j++)
3207 ir->opts.acc[i][j] = 0;
3211 nfrdim = str_nelem(frdim, MAXPTR, ptr1);
3212 nfreeze = str_nelem(freeze, MAXPTR, ptr2);
3213 if (nfrdim != DIM*nfreeze)
3215 gmx_fatal(FARGS, "Invalid Freezing input: %d groups and %d freeze values",
3218 do_numbering(natoms, groups, nfreeze, ptr2, grps, gnames, egcFREEZE,
3219 restnm, egrptpALL_GENREST, bVerbose, wi);
3220 nr = groups->grps[egcFREEZE].nr;
3221 ir->opts.ngfrz = nr;
3222 snew(ir->opts.nFreeze, nr);
3223 for (i = k = 0; (i < nfreeze); i++)
3225 for (j = 0; (j < DIM); j++, k++)
3227 ir->opts.nFreeze[i][j] = (gmx_strncasecmp(ptr1[k], "Y", 1) == 0);
3228 if (!ir->opts.nFreeze[i][j])
3230 if (gmx_strncasecmp(ptr1[k], "N", 1) != 0)
3232 sprintf(warnbuf, "Please use Y(ES) or N(O) for freezedim only "
3233 "(not %s)", ptr1[k]);
3234 warning(wi, warn_buf);
3239 for (; (i < nr); i++)
3241 for (j = 0; (j < DIM); j++)
3243 ir->opts.nFreeze[i][j] = 0;
3247 nenergy = str_nelem(energy, MAXPTR, ptr1);
3248 do_numbering(natoms, groups, nenergy, ptr1, grps, gnames, egcENER,
3249 restnm, egrptpALL_GENREST, bVerbose, wi);
3250 add_wall_energrps(groups, ir->nwall, symtab);
3251 ir->opts.ngener = groups->grps[egcENER].nr;
3252 nvcm = str_nelem(vcm, MAXPTR, ptr1);
3254 do_numbering(natoms, groups, nvcm, ptr1, grps, gnames, egcVCM,
3255 restnm, nvcm == 0 ? egrptpALL_GENREST : egrptpPART, bVerbose, wi);
3258 warning(wi, "Some atoms are not part of any center of mass motion removal group.\n"
3259 "This may lead to artifacts.\n"
3260 "In most cases one should use one group for the whole system.");
3263 /* Now we have filled the freeze struct, so we can calculate NRDF */
3264 calc_nrdf(mtop, ir, gnames);
3270 /* Must check per group! */
3271 for (i = 0; (i < ir->opts.ngtc); i++)
3273 ntot += ir->opts.nrdf[i];
3275 if (ntot != (DIM*natoms))
3277 fac = sqrt(ntot/(DIM*natoms));
3280 fprintf(stderr, "Scaling velocities by a factor of %.3f to account for constraints\n"
3281 "and removal of center of mass motion\n", fac);
3283 for (i = 0; (i < natoms); i++)
3285 svmul(fac, v[i], v[i]);
3290 nuser = str_nelem(user1, MAXPTR, ptr1);
3291 do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcUser1,
3292 restnm, egrptpALL_GENREST, bVerbose, wi);
3293 nuser = str_nelem(user2, MAXPTR, ptr1);
3294 do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcUser2,
3295 restnm, egrptpALL_GENREST, bVerbose, wi);
3296 nuser = str_nelem(xtc_grps, MAXPTR, ptr1);
3297 do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcXTC,
3298 restnm, egrptpONE, bVerbose, wi);
3299 nofg = str_nelem(orirefitgrp, MAXPTR, ptr1);
3300 do_numbering(natoms, groups, nofg, ptr1, grps, gnames, egcORFIT,
3301 restnm, egrptpALL_GENREST, bVerbose, wi);
3303 /* QMMM input processing */
3304 nQMg = str_nelem(QMMM, MAXPTR, ptr1);
3305 nQMmethod = str_nelem(QMmethod, MAXPTR, ptr2);
3306 nQMbasis = str_nelem(QMbasis, MAXPTR, ptr3);
3307 if ((nQMmethod != nQMg) || (nQMbasis != nQMg))
3309 gmx_fatal(FARGS, "Invalid QMMM input: %d groups %d basissets"
3310 " and %d methods\n", nQMg, nQMbasis, nQMmethod);
3312 /* group rest, if any, is always MM! */
3313 do_numbering(natoms, groups, nQMg, ptr1, grps, gnames, egcQMMM,
3314 restnm, egrptpALL_GENREST, bVerbose, wi);
3315 nr = nQMg; /*atoms->grps[egcQMMM].nr;*/
3316 ir->opts.ngQM = nQMg;
3317 snew(ir->opts.QMmethod, nr);
3318 snew(ir->opts.QMbasis, nr);
3319 for (i = 0; i < nr; i++)
3321 /* input consists of strings: RHF CASSCF PM3 .. These need to be
3322 * converted to the corresponding enum in names.c
3324 ir->opts.QMmethod[i] = search_QMstring(ptr2[i], eQMmethodNR,
3326 ir->opts.QMbasis[i] = search_QMstring(ptr3[i], eQMbasisNR,
3330 nQMmult = str_nelem(QMmult, MAXPTR, ptr1);
3331 nQMcharge = str_nelem(QMcharge, MAXPTR, ptr2);
3332 nbSH = str_nelem(bSH, MAXPTR, ptr3);
3333 snew(ir->opts.QMmult, nr);
3334 snew(ir->opts.QMcharge, nr);
3335 snew(ir->opts.bSH, nr);
3337 for (i = 0; i < nr; i++)
3339 ir->opts.QMmult[i] = strtol(ptr1[i], NULL, 10);
3340 ir->opts.QMcharge[i] = strtol(ptr2[i], NULL, 10);
3341 ir->opts.bSH[i] = (gmx_strncasecmp(ptr3[i], "Y", 1) == 0);
3344 nCASelec = str_nelem(CASelectrons, MAXPTR, ptr1);
3345 nCASorb = str_nelem(CASorbitals, MAXPTR, ptr2);
3346 snew(ir->opts.CASelectrons, nr);
3347 snew(ir->opts.CASorbitals, nr);
3348 for (i = 0; i < nr; i++)
3350 ir->opts.CASelectrons[i] = strtol(ptr1[i], NULL, 10);
3351 ir->opts.CASorbitals[i] = strtol(ptr2[i], NULL, 10);
3353 /* special optimization options */
3355 nbOPT = str_nelem(bOPT, MAXPTR, ptr1);
3356 nbTS = str_nelem(bTS, MAXPTR, ptr2);
3357 snew(ir->opts.bOPT, nr);
3358 snew(ir->opts.bTS, nr);
3359 for (i = 0; i < nr; i++)
3361 ir->opts.bOPT[i] = (gmx_strncasecmp(ptr1[i], "Y", 1) == 0);
3362 ir->opts.bTS[i] = (gmx_strncasecmp(ptr2[i], "Y", 1) == 0);
3364 nSAon = str_nelem(SAon, MAXPTR, ptr1);
3365 nSAoff = str_nelem(SAoff, MAXPTR, ptr2);
3366 nSAsteps = str_nelem(SAsteps, MAXPTR, ptr3);
3367 snew(ir->opts.SAon, nr);
3368 snew(ir->opts.SAoff, nr);
3369 snew(ir->opts.SAsteps, nr);
3371 for (i = 0; i < nr; i++)
3373 ir->opts.SAon[i] = strtod(ptr1[i], NULL);
3374 ir->opts.SAoff[i] = strtod(ptr2[i], NULL);
3375 ir->opts.SAsteps[i] = strtol(ptr3[i], NULL, 10);
3377 /* end of QMMM input */
3381 for (i = 0; (i < egcNR); i++)
3383 fprintf(stderr, "%-16s has %d element(s):", gtypes[i], groups->grps[i].nr);
3384 for (j = 0; (j < groups->grps[i].nr); j++)
3386 fprintf(stderr, " %s", *(groups->grpname[groups->grps[i].nm_ind[j]]));
3388 fprintf(stderr, "\n");
3392 nr = groups->grps[egcENER].nr;
3393 snew(ir->opts.egp_flags, nr*nr);
3395 bExcl = do_egp_flag(ir, groups, "energygrp-excl", egpexcl, EGP_EXCL);
3396 if (bExcl && ir->cutoff_scheme == ecutsVERLET)
3398 warning_error(wi, "Energy group exclusions are not (yet) implemented for the Verlet scheme");
3400 if (bExcl && EEL_FULL(ir->coulombtype))
3402 warning(wi, "Can not exclude the lattice Coulomb energy between energy groups");
3405 bTable = do_egp_flag(ir, groups, "energygrp-table", egptable, EGP_TABLE);
3406 if (bTable && !(ir->vdwtype == evdwUSER) &&
3407 !(ir->coulombtype == eelUSER) && !(ir->coulombtype == eelPMEUSER) &&
3408 !(ir->coulombtype == eelPMEUSERSWITCH))
3410 gmx_fatal(FARGS, "Can only have energy group pair tables in combination with user tables for VdW and/or Coulomb");
3413 decode_cos(efield_x, &(ir->ex[XX]));
3414 decode_cos(efield_xt, &(ir->et[XX]));
3415 decode_cos(efield_y, &(ir->ex[YY]));
3416 decode_cos(efield_yt, &(ir->et[YY]));
3417 decode_cos(efield_z, &(ir->ex[ZZ]));
3418 decode_cos(efield_zt, &(ir->et[ZZ]));
3422 do_adress_index(ir->adress, groups, gnames, &(ir->opts), wi);
3425 for (i = 0; (i < grps->nr); i++)
3437 static void check_disre(gmx_mtop_t *mtop)
3439 gmx_ffparams_t *ffparams;
3440 t_functype *functype;
3442 int i, ndouble, ftype;
3443 int label, old_label;
3445 if (gmx_mtop_ftype_count(mtop, F_DISRES) > 0)
3447 ffparams = &mtop->ffparams;
3448 functype = ffparams->functype;
3449 ip = ffparams->iparams;
3452 for (i = 0; i < ffparams->ntypes; i++)
3454 ftype = functype[i];
3455 if (ftype == F_DISRES)
3457 label = ip[i].disres.label;
3458 if (label == old_label)
3460 fprintf(stderr, "Distance restraint index %d occurs twice\n", label);
3468 gmx_fatal(FARGS, "Found %d double distance restraint indices,\n"
3469 "probably the parameters for multiple pairs in one restraint "
3470 "are not identical\n", ndouble);
3475 static gmx_bool absolute_reference(t_inputrec *ir, gmx_mtop_t *sys,
3476 gmx_bool posres_only,
3480 gmx_mtop_ilistloop_t iloop;
3490 for (d = 0; d < DIM; d++)
3492 AbsRef[d] = (d < ndof_com(ir) ? 0 : 1);
3494 /* Check for freeze groups */
3495 for (g = 0; g < ir->opts.ngfrz; g++)
3497 for (d = 0; d < DIM; d++)
3499 if (ir->opts.nFreeze[g][d] != 0)
3507 /* Check for position restraints */
3508 iloop = gmx_mtop_ilistloop_init(sys);
3509 while (gmx_mtop_ilistloop_next(iloop, &ilist, &nmol))
3512 (AbsRef[XX] == 0 || AbsRef[YY] == 0 || AbsRef[ZZ] == 0))
3514 for (i = 0; i < ilist[F_POSRES].nr; i += 2)
3516 pr = &sys->ffparams.iparams[ilist[F_POSRES].iatoms[i]];
3517 for (d = 0; d < DIM; d++)
3519 if (pr->posres.fcA[d] != 0)
3525 for (i = 0; i < ilist[F_FBPOSRES].nr; i += 2)
3527 /* Check for flat-bottom posres */
3528 pr = &sys->ffparams.iparams[ilist[F_FBPOSRES].iatoms[i]];
3529 if (pr->fbposres.k != 0)
3531 switch (pr->fbposres.geom)
3533 case efbposresSPHERE:
3534 AbsRef[XX] = AbsRef[YY] = AbsRef[ZZ] = 1;
3536 case efbposresCYLINDER:
3537 AbsRef[XX] = AbsRef[YY] = 1;
3539 case efbposresX: /* d=XX */
3540 case efbposresY: /* d=YY */
3541 case efbposresZ: /* d=ZZ */
3542 d = pr->fbposres.geom - efbposresX;
3546 gmx_fatal(FARGS, " Invalid geometry for flat-bottom position restraint.\n"
3547 "Expected nr between 1 and %d. Found %d\n", efbposresNR-1,
3555 return (AbsRef[XX] != 0 && AbsRef[YY] != 0 && AbsRef[ZZ] != 0);
3559 check_combination_rule_differences(const gmx_mtop_t *mtop, int state,
3560 gmx_bool *bC6ParametersWorkWithGeometricRules,
3561 gmx_bool *bC6ParametersWorkWithLBRules,
3562 gmx_bool *bLBRulesPossible)
3564 int ntypes, tpi, tpj, thisLBdiff, thisgeomdiff;
3567 #if (defined SIZEOF_LONG_LONG_INT) && (SIZEOF_LONG_LONG_INT >= 8)
3568 long long int npair, npair_ij;
3570 double npair, npair_ij;
3572 double geometricdiff, LBdiff;
3573 double c6i, c6j, c12i, c12j;
3574 double c6, c6_geometric, c6_LB;
3575 double sigmai, sigmaj, epsi, epsj;
3576 gmx_bool bCanDoLBRules, bCanDoGeometricRules;
3579 /* A tolerance of 1e-5 seems reasonable for (possibly hand-typed)
3580 * force-field floating point parameters.
3583 ptr = getenv("GMX_LJCOMB_TOL");
3588 sscanf(ptr, "%lf", &dbl);
3592 *bC6ParametersWorkWithLBRules = TRUE;
3593 *bC6ParametersWorkWithGeometricRules = TRUE;
3594 bCanDoLBRules = TRUE;
3595 bCanDoGeometricRules = TRUE;
3597 ntypes = mtop->ffparams.atnr;
3598 snew(typecount, ntypes);
3599 gmx_mtop_count_atomtypes(mtop, state, typecount);
3600 geometricdiff = LBdiff = 0.0;
3601 *bLBRulesPossible = TRUE;
3602 for (tpi = 0; tpi < ntypes; ++tpi)
3604 c6i = mtop->ffparams.iparams[(ntypes + 1) * tpi].lj.c6;
3605 c12i = mtop->ffparams.iparams[(ntypes + 1) * tpi].lj.c12;
3606 for (tpj = tpi; tpj < ntypes; ++tpj)
3608 c6j = mtop->ffparams.iparams[(ntypes + 1) * tpj].lj.c6;
3609 c12j = mtop->ffparams.iparams[(ntypes + 1) * tpj].lj.c12;
3610 c6 = mtop->ffparams.iparams[ntypes * tpi + tpj].lj.c6;
3611 c6_geometric = sqrt(c6i * c6j);
3612 if (!gmx_numzero(c6_geometric))
3614 if (!gmx_numzero(c12i) && !gmx_numzero(c12j))
3616 sigmai = pow(c12i / c6i, 1.0/6.0);
3617 sigmaj = pow(c12j / c6j, 1.0/6.0);
3618 epsi = c6i * c6i /(4.0 * c12i);
3619 epsj = c6j * c6j /(4.0 * c12j);
3620 c6_LB = 4.0 * pow(epsi * epsj, 1.0/2.0) * pow(0.5 * (sigmai + sigmaj), 6);
3624 *bLBRulesPossible = FALSE;
3625 c6_LB = c6_geometric;
3627 bCanDoLBRules = gmx_within_tol(c6_LB, c6, tol);
3630 if (FALSE == bCanDoLBRules)
3632 *bC6ParametersWorkWithLBRules = FALSE;
3635 bCanDoGeometricRules = gmx_within_tol(c6_geometric, c6, tol);
3637 if (FALSE == bCanDoGeometricRules)
3639 *bC6ParametersWorkWithGeometricRules = FALSE;
3647 check_combination_rules(const t_inputrec *ir, const gmx_mtop_t *mtop,
3651 gmx_bool bLBRulesPossible, bC6ParametersWorkWithGeometricRules, bC6ParametersWorkWithLBRules;
3653 check_combination_rule_differences(mtop, 0,
3654 &bC6ParametersWorkWithGeometricRules,
3655 &bC6ParametersWorkWithLBRules,
3657 if (ir->ljpme_combination_rule == eljpmeLB)
3659 if (FALSE == bC6ParametersWorkWithLBRules || FALSE == bLBRulesPossible)
3661 warning(wi, "You are using arithmetic-geometric combination rules "
3662 "in LJ-PME, but your non-bonded C6 parameters do not "
3663 "follow these rules.");
3668 if (FALSE == bC6ParametersWorkWithGeometricRules)
3670 if (ir->eDispCorr != edispcNO)
3672 warning_note(wi, "You are using geometric combination rules in "
3673 "LJ-PME, but your non-bonded C6 parameters do "
3674 "not follow these rules. "
3675 "If your force field uses Lorentz-Berthelot combination rules this "
3676 "will introduce small errors in the forces and energies in "
3677 "your simulations. Dispersion correction will correct total "
3678 "energy and/or pressure, but not forces or surface tensions. "
3679 "Please check the LJ-PME section in the manual "
3680 "before proceeding further.");
3684 warning_note(wi, "You are using geometric combination rules in "
3685 "LJ-PME, but your non-bonded C6 parameters do "
3686 "not follow these rules. "
3687 "If your force field uses Lorentz-Berthelot combination rules this "
3688 "will introduce small errors in the forces and energies in "
3689 "your simulations. Consider using dispersion correction "
3690 "for the total energy and pressure. "
3691 "Please check the LJ-PME section in the manual "
3692 "before proceeding further.");
3698 void triple_check(const char *mdparin, t_inputrec *ir, gmx_mtop_t *sys,
3702 int i, m, g, nmol, npct;
3703 gmx_bool bCharge, bAcc;
3704 real gdt_max, *mgrp, mt;
3706 gmx_mtop_atomloop_block_t aloopb;
3707 gmx_mtop_atomloop_all_t aloop;
3710 char warn_buf[STRLEN];
3712 set_warning_line(wi, mdparin, -1);
3714 if (EI_DYNAMICS(ir->eI) && !EI_SD(ir->eI) && ir->eI != eiBD &&
3715 ir->comm_mode == ecmNO &&
3716 !(absolute_reference(ir, sys, FALSE, AbsRef) || ir->nsteps <= 10))
3718 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");
3721 /* Check for pressure coupling with absolute position restraints */
3722 if (ir->epc != epcNO && ir->refcoord_scaling == erscNO)
3724 absolute_reference(ir, sys, TRUE, AbsRef);
3726 for (m = 0; m < DIM; m++)
3728 if (AbsRef[m] && norm2(ir->compress[m]) > 0)
3730 warning(wi, "You are using pressure coupling with absolute position restraints, this will give artifacts. Use the refcoord_scaling option.");
3738 aloopb = gmx_mtop_atomloop_block_init(sys);
3739 while (gmx_mtop_atomloop_block_next(aloopb, &atom, &nmol))
3741 if (atom->q != 0 || atom->qB != 0)
3749 if (EEL_FULL(ir->coulombtype))
3752 "You are using full electrostatics treatment %s for a system without charges.\n"
3753 "This costs a lot of performance for just processing zeros, consider using %s instead.\n",
3754 EELTYPE(ir->coulombtype), EELTYPE(eelCUT));
3755 warning(wi, err_buf);
3760 if (ir->coulombtype == eelCUT && ir->rcoulomb > 0 && !ir->implicit_solvent)
3763 "You are using a plain Coulomb cut-off, which might produce artifacts.\n"
3764 "You might want to consider using %s electrostatics.\n",
3766 warning_note(wi, err_buf);
3770 /* Check if combination rules used in LJ-PME are the same as in the force field */
3771 if (EVDW_PME(ir->vdwtype))
3773 check_combination_rules(ir, sys, wi);
3776 /* Generalized reaction field */
3777 if (ir->opts.ngtc == 0)
3779 sprintf(err_buf, "No temperature coupling while using coulombtype %s",
3781 CHECK(ir->coulombtype == eelGRF);
3785 sprintf(err_buf, "When using coulombtype = %s"
3786 " ref-t for temperature coupling should be > 0",
3788 CHECK((ir->coulombtype == eelGRF) && (ir->opts.ref_t[0] <= 0));
3791 if (ir->eI == eiSD1 &&
3792 (gmx_mtop_ftype_count(sys, F_CONSTR) > 0 ||
3793 gmx_mtop_ftype_count(sys, F_SETTLE) > 0))
3795 sprintf(warn_buf, "With constraints integrator %s is less accurate, consider using %s instead", ei_names[ir->eI], ei_names[eiSD2]);
3796 warning_note(wi, warn_buf);
3800 for (i = 0; (i < sys->groups.grps[egcACC].nr); i++)
3802 for (m = 0; (m < DIM); m++)
3804 if (fabs(ir->opts.acc[i][m]) > 1e-6)
3813 snew(mgrp, sys->groups.grps[egcACC].nr);
3814 aloop = gmx_mtop_atomloop_all_init(sys);
3815 while (gmx_mtop_atomloop_all_next(aloop, &i, &atom))
3817 mgrp[ggrpnr(&sys->groups, egcACC, i)] += atom->m;
3820 for (i = 0; (i < sys->groups.grps[egcACC].nr); i++)
3822 for (m = 0; (m < DIM); m++)
3824 acc[m] += ir->opts.acc[i][m]*mgrp[i];
3828 for (m = 0; (m < DIM); m++)
3830 if (fabs(acc[m]) > 1e-6)
3832 const char *dim[DIM] = { "X", "Y", "Z" };
3834 "Net Acceleration in %s direction, will %s be corrected\n",
3835 dim[m], ir->nstcomm != 0 ? "" : "not");
3836 if (ir->nstcomm != 0 && m < ndof_com(ir))
3839 for (i = 0; (i < sys->groups.grps[egcACC].nr); i++)
3841 ir->opts.acc[i][m] -= acc[m];
3849 if (ir->efep != efepNO && ir->fepvals->sc_alpha != 0 &&
3850 !gmx_within_tol(sys->ffparams.reppow, 12.0, 10*GMX_DOUBLE_EPS))
3852 gmx_fatal(FARGS, "Soft-core interactions are only supported with VdW repulsion power 12");
3855 if (ir->ePull != epullNO)
3857 if (ir->pull->grp[0].nat == 0)
3859 absolute_reference(ir, sys, FALSE, AbsRef);
3860 for (m = 0; m < DIM; m++)
3862 if (ir->pull->dim[m] && !AbsRef[m])
3864 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.");
3870 if (ir->pull->eGeom == epullgDIRPBC)
3872 for (i = 0; i < 3; i++)
3874 for (m = 0; m <= i; m++)
3876 if ((ir->epc != epcNO && ir->compress[i][m] != 0) ||
3877 ir->deform[i][m] != 0)
3879 for (g = 1; g < ir->pull->ngrp; g++)
3881 if (ir->pull->grp[g].vec[m] != 0)
3883 gmx_fatal(FARGS, "Can not have dynamic box while using pull geometry '%s' (dim %c)", EPULLGEOM(ir->pull->eGeom), 'x'+m);
3895 void double_check(t_inputrec *ir, matrix box, gmx_bool bConstr, warninp_t wi)
3899 char warn_buf[STRLEN];
3902 ptr = check_box(ir->ePBC, box);
3905 warning_error(wi, ptr);
3908 if (bConstr && ir->eConstrAlg == econtSHAKE)
3910 if (ir->shake_tol <= 0.0)
3912 sprintf(warn_buf, "ERROR: shake-tol must be > 0 instead of %g\n",
3914 warning_error(wi, warn_buf);
3917 if (IR_TWINRANGE(*ir) && ir->nstlist > 1)
3919 sprintf(warn_buf, "With twin-range cut-off's and SHAKE the virial and the pressure are incorrect.");
3920 if (ir->epc == epcNO)
3922 warning(wi, warn_buf);
3926 warning_error(wi, warn_buf);
3931 if ( (ir->eConstrAlg == econtLINCS) && bConstr)
3933 /* If we have Lincs constraints: */
3934 if (ir->eI == eiMD && ir->etc == etcNO &&
3935 ir->eConstrAlg == econtLINCS && ir->nLincsIter == 1)
3937 sprintf(warn_buf, "For energy conservation with LINCS, lincs_iter should be 2 or larger.\n");
3938 warning_note(wi, warn_buf);
3941 if ((ir->eI == eiCG || ir->eI == eiLBFGS) && (ir->nProjOrder < 8))
3943 sprintf(warn_buf, "For accurate %s with LINCS constraints, lincs-order should be 8 or more.", ei_names[ir->eI]);
3944 warning_note(wi, warn_buf);
3946 if (ir->epc == epcMTTK)
3948 warning_error(wi, "MTTK not compatible with lincs -- use shake instead.");
3952 if (ir->LincsWarnAngle > 90.0)
3954 sprintf(warn_buf, "lincs-warnangle can not be larger than 90 degrees, setting it to 90.\n");
3955 warning(wi, warn_buf);
3956 ir->LincsWarnAngle = 90.0;
3959 if (ir->ePBC != epbcNONE)
3961 if (ir->nstlist == 0)
3963 warning(wi, "With nstlist=0 atoms are only put into the box at step 0, therefore drifting atoms might cause the simulation to crash.");
3965 bTWIN = (ir->rlistlong > ir->rlist);
3966 if (ir->ns_type == ensGRID)
3968 if (sqr(ir->rlistlong) >= max_cutoff2(ir->ePBC, box))
3970 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",
3971 bTWIN ? (ir->rcoulomb == ir->rlistlong ? "rcoulomb" : "rvdw") : "rlist");
3972 warning_error(wi, warn_buf);
3977 min_size = min(box[XX][XX], min(box[YY][YY], box[ZZ][ZZ]));
3978 if (2*ir->rlistlong >= min_size)
3980 sprintf(warn_buf, "ERROR: One of the box lengths is smaller than twice the cut-off length. Increase the box size or decrease rlist.");
3981 warning_error(wi, warn_buf);
3984 fprintf(stderr, "Grid search might allow larger cut-off's than simple search with triclinic boxes.");
3991 void check_chargegroup_radii(const gmx_mtop_t *mtop, const t_inputrec *ir,
3995 real rvdw1, rvdw2, rcoul1, rcoul2;
3996 char warn_buf[STRLEN];
3998 calc_chargegroup_radii(mtop, x, &rvdw1, &rvdw2, &rcoul1, &rcoul2);
4002 printf("Largest charge group radii for Van der Waals: %5.3f, %5.3f nm\n",
4007 printf("Largest charge group radii for Coulomb: %5.3f, %5.3f nm\n",
4013 if (rvdw1 + rvdw2 > ir->rlist ||
4014 rcoul1 + rcoul2 > ir->rlist)
4017 "The sum of the two largest charge group radii (%f) "
4018 "is larger than rlist (%f)\n",
4019 max(rvdw1+rvdw2, rcoul1+rcoul2), ir->rlist);
4020 warning(wi, warn_buf);
4024 /* Here we do not use the zero at cut-off macro,
4025 * since user defined interactions might purposely
4026 * not be zero at the cut-off.
4028 if ((EVDW_IS_ZERO_AT_CUTOFF(ir->vdwtype) ||
4029 ir->vdw_modifier != eintmodNONE) &&
4030 rvdw1 + rvdw2 > ir->rlistlong - ir->rvdw)
4032 sprintf(warn_buf, "The sum of the two largest charge group "
4033 "radii (%f) is larger than %s (%f) - rvdw (%f).\n"
4034 "With exact cut-offs, better performance can be "
4035 "obtained with cutoff-scheme = %s, because it "
4036 "does not use charge groups at all.",
4038 ir->rlistlong > ir->rlist ? "rlistlong" : "rlist",
4039 ir->rlistlong, ir->rvdw,
4040 ecutscheme_names[ecutsVERLET]);
4043 warning(wi, warn_buf);
4047 warning_note(wi, warn_buf);
4050 if ((EEL_IS_ZERO_AT_CUTOFF(ir->coulombtype) ||
4051 ir->coulomb_modifier != eintmodNONE) &&
4052 rcoul1 + rcoul2 > ir->rlistlong - ir->rcoulomb)
4054 sprintf(warn_buf, "The sum of the two largest charge group radii (%f) is larger than %s (%f) - rcoulomb (%f).\n"
4055 "With exact cut-offs, better performance can be obtained with cutoff-scheme = %s, because it does not use charge groups at all.",
4057 ir->rlistlong > ir->rlist ? "rlistlong" : "rlist",
4058 ir->rlistlong, ir->rcoulomb,
4059 ecutscheme_names[ecutsVERLET]);
4062 warning(wi, warn_buf);
4066 warning_note(wi, warn_buf);