2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
49 #include "gmx_fatal.h"
62 #include "mtop_util.h"
63 #include "chargegroup.h"
65 #include "calc_verletbuf.h"
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 typedef struct t_inputrec_strings
79 char tcgrps[STRLEN], tau_t[STRLEN], ref_t[STRLEN],
80 acc[STRLEN], accgrps[STRLEN], freeze[STRLEN], frdim[STRLEN],
81 energy[STRLEN], user1[STRLEN], user2[STRLEN], vcm[STRLEN], x_compressed_groups[STRLEN],
82 couple_moltype[STRLEN], orirefitgrp[STRLEN], egptable[STRLEN], egpexcl[STRLEN],
83 wall_atomtype[STRLEN], wall_density[STRLEN], deform[STRLEN], QMMM[STRLEN];
84 char fep_lambda[efptNR][STRLEN];
85 char lambda_weights[STRLEN];
88 char anneal[STRLEN], anneal_npoints[STRLEN],
89 anneal_time[STRLEN], anneal_temp[STRLEN];
90 char QMmethod[STRLEN], QMbasis[STRLEN], QMcharge[STRLEN], QMmult[STRLEN],
91 bSH[STRLEN], CASorbitals[STRLEN], CASelectrons[STRLEN], SAon[STRLEN],
92 SAoff[STRLEN], SAsteps[STRLEN], bTS[STRLEN], bOPT[STRLEN];
93 char efield_x[STRLEN], efield_xt[STRLEN], efield_y[STRLEN],
94 efield_yt[STRLEN], efield_z[STRLEN], efield_zt[STRLEN];
96 } gmx_inputrec_strings;
98 static gmx_inputrec_strings *is = NULL;
100 void init_inputrec_strings()
104 gmx_incons("Attempted to call init_inputrec_strings before calling done_inputrec_strings. Only one inputrec (i.e. .mdp file) can be parsed at a time.");
109 void done_inputrec_strings()
115 static char swapgrp[STRLEN], splitgrp0[STRLEN], splitgrp1[STRLEN], solgrp[STRLEN];
118 egrptpALL, /* All particles have to be a member of a group. */
119 egrptpALL_GENREST, /* A rest group with name is generated for particles *
120 * that are not part of any group. */
121 egrptpPART, /* As egrptpALL_GENREST, but no name is generated *
122 * for the rest group. */
123 egrptpONE /* Merge all selected groups into one group, *
124 * make a rest group for the remaining particles. */
127 static const char *constraints[eshNR+1] = {
128 "none", "h-bonds", "all-bonds", "h-angles", "all-angles", NULL
131 static const char *couple_lam[ecouplamNR+1] = {
132 "vdw-q", "vdw", "q", "none", NULL
135 void init_ir(t_inputrec *ir, t_gromppopts *opts)
137 snew(opts->include, STRLEN);
138 snew(opts->define, STRLEN);
139 snew(ir->fepvals, 1);
140 snew(ir->expandedvals, 1);
141 snew(ir->simtempvals, 1);
144 static void GetSimTemps(int ntemps, t_simtemp *simtemp, double *temperature_lambdas)
149 for (i = 0; i < ntemps; i++)
151 /* simple linear scaling -- allows more control */
152 if (simtemp->eSimTempScale == esimtempLINEAR)
154 simtemp->temperatures[i] = simtemp->simtemp_low + (simtemp->simtemp_high-simtemp->simtemp_low)*temperature_lambdas[i];
156 else if (simtemp->eSimTempScale == esimtempGEOMETRIC) /* should give roughly equal acceptance for constant heat capacity . . . */
158 simtemp->temperatures[i] = simtemp->simtemp_low * pow(simtemp->simtemp_high/simtemp->simtemp_low, (1.0*i)/(ntemps-1));
160 else if (simtemp->eSimTempScale == esimtempEXPONENTIAL)
162 simtemp->temperatures[i] = simtemp->simtemp_low + (simtemp->simtemp_high-simtemp->simtemp_low)*((exp(temperature_lambdas[i])-1)/(exp(1.0)-1));
167 sprintf(errorstr, "eSimTempScale=%d not defined", simtemp->eSimTempScale);
168 gmx_fatal(FARGS, errorstr);
175 static void _low_check(gmx_bool b, char *s, warninp_t wi)
179 warning_error(wi, s);
183 static void check_nst(const char *desc_nst, int nst,
184 const char *desc_p, int *p,
189 if (*p > 0 && *p % nst != 0)
191 /* Round up to the next multiple of nst */
192 *p = ((*p)/nst + 1)*nst;
193 sprintf(buf, "%s should be a multiple of %s, changing %s to %d\n",
194 desc_p, desc_nst, desc_p, *p);
199 static gmx_bool ir_NVE(const t_inputrec *ir)
201 return ((ir->eI == eiMD || EI_VV(ir->eI)) && ir->etc == etcNO);
204 static int lcd(int n1, int n2)
209 for (i = 2; (i <= n1 && i <= n2); i++)
211 if (n1 % i == 0 && n2 % i == 0)
220 static void process_interaction_modifier(const t_inputrec *ir, int *eintmod)
222 if (*eintmod == eintmodPOTSHIFT_VERLET)
224 if (ir->cutoff_scheme == ecutsVERLET)
226 *eintmod = eintmodPOTSHIFT;
230 *eintmod = eintmodNONE;
235 void check_ir(const char *mdparin, t_inputrec *ir, t_gromppopts *opts,
237 /* Check internal consistency */
239 /* Strange macro: first one fills the err_buf, and then one can check
240 * the condition, which will print the message and increase the error
243 #define CHECK(b) _low_check(b, err_buf, wi)
244 char err_buf[256], warn_buf[STRLEN];
250 t_lambda *fep = ir->fepvals;
251 t_expanded *expand = ir->expandedvals;
253 set_warning_line(wi, mdparin, -1);
255 /* BASIC CUT-OFF STUFF */
256 if (ir->rcoulomb < 0)
258 warning_error(wi, "rcoulomb should be >= 0");
262 warning_error(wi, "rvdw should be >= 0");
265 !(ir->cutoff_scheme == ecutsVERLET && ir->verletbuf_tol > 0))
267 warning_error(wi, "rlist should be >= 0");
270 process_interaction_modifier(ir, &ir->coulomb_modifier);
271 process_interaction_modifier(ir, &ir->vdw_modifier);
273 if (ir->cutoff_scheme == ecutsGROUP)
276 "The group cutoff scheme is deprecated in Gromacs 5.0 and will be removed in a future "
277 "release when all interaction forms are supported for the verlet scheme. The verlet "
278 "scheme already scales better, and it is compatible with GPUs and other accelerators.");
280 /* BASIC CUT-OFF STUFF */
281 if (ir->rlist == 0 ||
282 !((ir_coulomb_might_be_zero_at_cutoff(ir) && ir->rcoulomb > ir->rlist) ||
283 (ir_vdw_might_be_zero_at_cutoff(ir) && ir->rvdw > ir->rlist)))
285 /* No switched potential and/or no twin-range:
286 * we can set the long-range cut-off to the maximum of the other cut-offs.
288 ir->rlistlong = max_cutoff(ir->rlist, max_cutoff(ir->rvdw, ir->rcoulomb));
290 else if (ir->rlistlong < 0)
292 ir->rlistlong = max_cutoff(ir->rlist, max_cutoff(ir->rvdw, ir->rcoulomb));
293 sprintf(warn_buf, "rlistlong was not set, setting it to %g (no buffer)",
295 warning(wi, warn_buf);
297 if (ir->rlistlong == 0 && ir->ePBC != epbcNONE)
299 warning_error(wi, "Can not have an infinite cut-off with PBC");
301 if (ir->rlistlong > 0 && (ir->rlist == 0 || ir->rlistlong < ir->rlist))
303 warning_error(wi, "rlistlong can not be shorter than rlist");
305 if (IR_TWINRANGE(*ir) && ir->nstlist <= 0)
307 warning_error(wi, "Can not have nstlist<=0 with twin-range interactions");
311 if (ir->rlistlong == ir->rlist)
315 else if (ir->rlistlong > ir->rlist && ir->nstcalclr == 0)
317 warning_error(wi, "With different cutoffs for electrostatics and VdW, nstcalclr must be -1 or a positive number");
320 if (ir->cutoff_scheme == ecutsVERLET)
324 /* Normal Verlet type neighbor-list, currently only limited feature support */
325 if (inputrec2nboundeddim(ir) < 3)
327 warning_error(wi, "With Verlet lists only full pbc or pbc=xy with walls is supported");
329 if (ir->rcoulomb != ir->rvdw)
331 warning_error(wi, "With Verlet lists rcoulomb!=rvdw is not supported");
333 if (ir->vdwtype == evdwSHIFT || ir->vdwtype == evdwSWITCH)
335 if (ir->vdw_modifier == eintmodNONE ||
336 ir->vdw_modifier == eintmodPOTSHIFT)
338 ir->vdw_modifier = (ir->vdwtype == evdwSHIFT ? eintmodFORCESWITCH : eintmodPOTSWITCH);
340 sprintf(warn_buf, "Replacing vdwtype=%s by the equivalent combination of vdwtype=%s and vdw_modifier=%s", evdw_names[ir->vdwtype], evdw_names[evdwCUT], eintmod_names[ir->vdw_modifier]);
341 warning_note(wi, warn_buf);
343 ir->vdwtype = evdwCUT;
347 sprintf(warn_buf, "Unsupported combination of vdwtype=%s and vdw_modifier=%s", evdw_names[ir->vdwtype], eintmod_names[ir->vdw_modifier]);
348 warning_error(wi, warn_buf);
352 if (!(ir->vdwtype == evdwCUT || ir->vdwtype == evdwPME))
354 warning_error(wi, "With Verlet lists only cut-off and PME LJ interactions are supported");
356 if (!(ir->coulombtype == eelCUT ||
357 (EEL_RF(ir->coulombtype) && ir->coulombtype != eelRF_NEC) ||
358 EEL_PME(ir->coulombtype) || ir->coulombtype == eelEWALD))
360 warning_error(wi, "With Verlet lists only cut-off, reaction-field, PME and Ewald electrostatics are supported");
362 if (!(ir->coulomb_modifier == eintmodNONE ||
363 ir->coulomb_modifier == eintmodPOTSHIFT))
365 sprintf(warn_buf, "coulomb_modifier=%s is not supported with the Verlet cut-off scheme", eintmod_names[ir->coulomb_modifier]);
366 warning_error(wi, warn_buf);
369 if (ir->nstlist <= 0)
371 warning_error(wi, "With Verlet lists nstlist should be larger than 0");
374 if (ir->nstlist < 10)
376 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.");
379 rc_max = max(ir->rvdw, ir->rcoulomb);
381 if (ir->verletbuf_tol <= 0)
383 if (ir->verletbuf_tol == 0)
385 warning_error(wi, "Can not have Verlet buffer tolerance of exactly 0");
388 if (ir->rlist < rc_max)
390 warning_error(wi, "With verlet lists rlist can not be smaller than rvdw or rcoulomb");
393 if (ir->rlist == rc_max && ir->nstlist > 1)
395 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.");
400 if (ir->rlist > rc_max)
402 warning_note(wi, "You have set rlist larger than the interaction cut-off, but you also have verlet-buffer-tolerance > 0. Will set rlist using verlet-buffer-tolerance.");
405 if (ir->nstlist == 1)
407 /* No buffer required */
412 if (EI_DYNAMICS(ir->eI))
414 if (inputrec2nboundeddim(ir) < 3)
416 warning_error(wi, "The box volume is required for calculating rlist from the energy drift with verlet-buffer-tolerance > 0. You are using at least one unbounded dimension, so no volume can be computed. Either use a finite box, or set rlist yourself together with verlet-buffer-tolerance = -1.");
418 /* Set rlist temporarily so we can continue processing */
423 /* Set the buffer to 5% of the cut-off */
424 ir->rlist = (1.0 + verlet_buffer_ratio_nodynamics)*rc_max;
429 /* No twin-range calculations with Verlet lists */
430 ir->rlistlong = ir->rlist;
433 if (ir->nstcalclr == -1)
435 /* if rlist=rlistlong, this will later be changed to nstcalclr=0 */
436 ir->nstcalclr = ir->nstlist;
438 else if (ir->nstcalclr > 0)
440 if (ir->nstlist > 0 && (ir->nstlist % ir->nstcalclr != 0))
442 warning_error(wi, "nstlist must be evenly divisible by nstcalclr. Use nstcalclr = -1 to automatically follow nstlist");
445 else if (ir->nstcalclr < -1)
447 warning_error(wi, "nstcalclr must be a positive number (divisor of nstcalclr), or -1 to follow nstlist.");
450 if (EEL_PME(ir->coulombtype) && ir->rcoulomb > ir->rvdw && ir->nstcalclr > 1)
452 warning_error(wi, "When used with PME, the long-range component of twin-range interactions must be updated every step (nstcalclr)");
455 /* GENERAL INTEGRATOR STUFF */
456 if (!(ir->eI == eiMD || EI_VV(ir->eI)))
460 if (ir->eI == eiVVAK)
462 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]);
463 warning_note(wi, warn_buf);
465 if (!EI_DYNAMICS(ir->eI))
469 if (EI_DYNAMICS(ir->eI))
471 if (ir->nstcalcenergy < 0)
473 ir->nstcalcenergy = ir_optimal_nstcalcenergy(ir);
474 if (ir->nstenergy != 0 && ir->nstenergy < ir->nstcalcenergy)
476 /* nstcalcenergy larger than nstener does not make sense.
477 * We ideally want nstcalcenergy=nstener.
481 ir->nstcalcenergy = lcd(ir->nstenergy, ir->nstlist);
485 ir->nstcalcenergy = ir->nstenergy;
489 else if ( (ir->nstenergy > 0 && ir->nstcalcenergy > ir->nstenergy) ||
490 (ir->efep != efepNO && ir->fepvals->nstdhdl > 0 &&
491 (ir->nstcalcenergy > ir->fepvals->nstdhdl) ) )
494 const char *nsten = "nstenergy";
495 const char *nstdh = "nstdhdl";
496 const char *min_name = nsten;
497 int min_nst = ir->nstenergy;
499 /* find the smallest of ( nstenergy, nstdhdl ) */
500 if (ir->efep != efepNO && ir->fepvals->nstdhdl > 0 &&
501 (ir->nstenergy == 0 || ir->fepvals->nstdhdl < ir->nstenergy))
503 min_nst = ir->fepvals->nstdhdl;
506 /* If the user sets nstenergy small, we should respect that */
508 "Setting nstcalcenergy (%d) equal to %s (%d)",
509 ir->nstcalcenergy, min_name, min_nst);
510 warning_note(wi, warn_buf);
511 ir->nstcalcenergy = min_nst;
514 if (ir->epc != epcNO)
516 if (ir->nstpcouple < 0)
518 ir->nstpcouple = ir_optimal_nstpcouple(ir);
521 if (IR_TWINRANGE(*ir))
523 check_nst("nstlist", ir->nstlist,
524 "nstcalcenergy", &ir->nstcalcenergy, wi);
525 if (ir->epc != epcNO)
527 check_nst("nstlist", ir->nstlist,
528 "nstpcouple", &ir->nstpcouple, wi);
532 if (ir->nstcalcenergy > 0)
534 if (ir->efep != efepNO)
536 /* nstdhdl should be a multiple of nstcalcenergy */
537 check_nst("nstcalcenergy", ir->nstcalcenergy,
538 "nstdhdl", &ir->fepvals->nstdhdl, wi);
539 /* nstexpanded should be a multiple of nstcalcenergy */
540 check_nst("nstcalcenergy", ir->nstcalcenergy,
541 "nstexpanded", &ir->expandedvals->nstexpanded, wi);
543 /* for storing exact averages nstenergy should be
544 * a multiple of nstcalcenergy
546 check_nst("nstcalcenergy", ir->nstcalcenergy,
547 "nstenergy", &ir->nstenergy, wi);
552 if ((EI_SD(ir->eI) || ir->eI == eiBD) &&
553 ir->bContinuation && ir->ld_seed != -1)
555 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)");
561 sprintf(err_buf, "TPI only works with pbc = %s", epbc_names[epbcXYZ]);
562 CHECK(ir->ePBC != epbcXYZ);
563 sprintf(err_buf, "TPI only works with ns = %s", ens_names[ensGRID]);
564 CHECK(ir->ns_type != ensGRID);
565 sprintf(err_buf, "with TPI nstlist should be larger than zero");
566 CHECK(ir->nstlist <= 0);
567 sprintf(err_buf, "TPI does not work with full electrostatics other than PME");
568 CHECK(EEL_FULL(ir->coulombtype) && !EEL_PME(ir->coulombtype));
572 if ( (opts->nshake > 0) && (opts->bMorse) )
575 "Using morse bond-potentials while constraining bonds is useless");
576 warning(wi, warn_buf);
579 if ((EI_SD(ir->eI) || ir->eI == eiBD) &&
580 ir->bContinuation && ir->ld_seed != -1)
582 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)");
584 /* verify simulated tempering options */
588 gmx_bool bAllTempZero = TRUE;
589 for (i = 0; i < fep->n_lambda; i++)
591 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]);
592 CHECK((fep->all_lambda[efptTEMPERATURE][i] < 0) || (fep->all_lambda[efptTEMPERATURE][i] > 1));
593 if (fep->all_lambda[efptTEMPERATURE][i] > 0)
595 bAllTempZero = FALSE;
598 sprintf(err_buf, "if simulated tempering is on, temperature-lambdas may not be all zero");
599 CHECK(bAllTempZero == TRUE);
601 sprintf(err_buf, "Simulated tempering is currently only compatible with md-vv");
602 CHECK(ir->eI != eiVV);
604 /* check compatability of the temperature coupling with simulated tempering */
606 if (ir->etc == etcNOSEHOOVER)
608 sprintf(warn_buf, "Nose-Hoover based temperature control such as [%s] my not be entirelyconsistent with simulated tempering", etcoupl_names[ir->etc]);
609 warning_note(wi, warn_buf);
612 /* check that the temperatures make sense */
614 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);
615 CHECK(ir->simtempvals->simtemp_high <= ir->simtempvals->simtemp_low);
617 sprintf(err_buf, "Higher simulated tempering temperature (%g) must be >= zero", ir->simtempvals->simtemp_high);
618 CHECK(ir->simtempvals->simtemp_high <= 0);
620 sprintf(err_buf, "Lower simulated tempering temperature (%g) must be >= zero", ir->simtempvals->simtemp_low);
621 CHECK(ir->simtempvals->simtemp_low <= 0);
624 /* verify free energy options */
626 if (ir->efep != efepNO)
629 sprintf(err_buf, "The soft-core power is %d and can only be 1 or 2",
631 CHECK(fep->sc_alpha != 0 && fep->sc_power != 1 && fep->sc_power != 2);
633 sprintf(err_buf, "The soft-core sc-r-power is %d and can only be 6 or 48",
634 (int)fep->sc_r_power);
635 CHECK(fep->sc_alpha != 0 && fep->sc_r_power != 6.0 && fep->sc_r_power != 48.0);
637 sprintf(err_buf, "Can't use postive delta-lambda (%g) if initial state/lambda does not start at zero", fep->delta_lambda);
638 CHECK(fep->delta_lambda > 0 && ((fep->init_fep_state > 0) || (fep->init_lambda > 0)));
640 sprintf(err_buf, "Can't use postive delta-lambda (%g) with expanded ensemble simulations", fep->delta_lambda);
641 CHECK(fep->delta_lambda > 0 && (ir->efep == efepEXPANDED));
643 sprintf(err_buf, "Can only use expanded ensemble with md-vv for now; should be supported for other integrators in 5.0");
644 CHECK(!(EI_VV(ir->eI)) && (ir->efep == efepEXPANDED));
646 sprintf(err_buf, "Free-energy not implemented for Ewald");
647 CHECK(ir->coulombtype == eelEWALD);
649 /* check validty of lambda inputs */
650 if (fep->n_lambda == 0)
652 /* Clear output in case of no states:*/
653 sprintf(err_buf, "init-lambda-state set to %d: no lambda states are defined.", fep->init_fep_state);
654 CHECK((fep->init_fep_state >= 0) && (fep->n_lambda == 0));
658 sprintf(err_buf, "initial thermodynamic state %d does not exist, only goes to %d", fep->init_fep_state, fep->n_lambda-1);
659 CHECK((fep->init_fep_state >= fep->n_lambda));
662 sprintf(err_buf, "Lambda state must be set, either with init-lambda-state or with init-lambda");
663 CHECK((fep->init_fep_state < 0) && (fep->init_lambda < 0));
665 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",
666 fep->init_lambda, fep->init_fep_state);
667 CHECK((fep->init_fep_state >= 0) && (fep->init_lambda >= 0));
671 if ((fep->init_lambda >= 0) && (fep->delta_lambda == 0))
675 for (i = 0; i < efptNR; i++)
677 if (fep->separate_dvdl[i])
682 if (n_lambda_terms > 1)
684 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.");
685 warning(wi, warn_buf);
688 if (n_lambda_terms < 2 && fep->n_lambda > 0)
691 "init-lambda is deprecated for setting lambda state (except for slow growth). Use init-lambda-state instead.");
695 for (j = 0; j < efptNR; j++)
697 for (i = 0; i < fep->n_lambda; i++)
699 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]);
700 CHECK((fep->all_lambda[j][i] < 0) || (fep->all_lambda[j][i] > 1));
704 if ((fep->sc_alpha > 0) && (!fep->bScCoul))
706 for (i = 0; i < fep->n_lambda; i++)
708 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],
709 fep->all_lambda[efptCOUL][i]);
710 CHECK((fep->sc_alpha > 0) &&
711 (((fep->all_lambda[efptCOUL][i] > 0.0) &&
712 (fep->all_lambda[efptCOUL][i] < 1.0)) &&
713 ((fep->all_lambda[efptVDW][i] > 0.0) &&
714 (fep->all_lambda[efptVDW][i] < 1.0))));
718 if ((fep->bScCoul) && (EEL_PME(ir->coulombtype)))
720 real sigma, lambda, r_sc;
723 /* Maximum estimate for A and B charges equal with lambda power 1 */
725 r_sc = pow(lambda*fep->sc_alpha*pow(sigma/ir->rcoulomb, fep->sc_r_power) + 1.0, 1.0/fep->sc_r_power);
726 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.",
728 sigma, lambda, r_sc - 1.0, ir->ewald_rtol);
729 warning_note(wi, warn_buf);
732 /* Free Energy Checks -- In an ideal world, slow growth and FEP would
733 be treated differently, but that's the next step */
735 for (i = 0; i < efptNR; i++)
737 for (j = 0; j < fep->n_lambda; j++)
739 sprintf(err_buf, "%s[%d] must be between 0 and 1", efpt_names[i], j);
740 CHECK((fep->all_lambda[i][j] < 0) || (fep->all_lambda[i][j] > 1));
745 if ((ir->bSimTemp) || (ir->efep == efepEXPANDED))
748 expand = ir->expandedvals;
750 /* checking equilibration of weights inputs for validity */
752 sprintf(err_buf, "weight-equil-number-all-lambda (%d) is ignored if lmc-weights-equil is not equal to %s",
753 expand->equil_n_at_lam, elmceq_names[elmceqNUMATLAM]);
754 CHECK((expand->equil_n_at_lam > 0) && (expand->elmceq != elmceqNUMATLAM));
756 sprintf(err_buf, "weight-equil-number-samples (%d) is ignored if lmc-weights-equil is not equal to %s",
757 expand->equil_samples, elmceq_names[elmceqSAMPLES]);
758 CHECK((expand->equil_samples > 0) && (expand->elmceq != elmceqSAMPLES));
760 sprintf(err_buf, "weight-equil-number-steps (%d) is ignored if lmc-weights-equil is not equal to %s",
761 expand->equil_steps, elmceq_names[elmceqSTEPS]);
762 CHECK((expand->equil_steps > 0) && (expand->elmceq != elmceqSTEPS));
764 sprintf(err_buf, "weight-equil-wl-delta (%d) is ignored if lmc-weights-equil is not equal to %s",
765 expand->equil_samples, elmceq_names[elmceqWLDELTA]);
766 CHECK((expand->equil_wl_delta > 0) && (expand->elmceq != elmceqWLDELTA));
768 sprintf(err_buf, "weight-equil-count-ratio (%f) is ignored if lmc-weights-equil is not equal to %s",
769 expand->equil_ratio, elmceq_names[elmceqRATIO]);
770 CHECK((expand->equil_ratio > 0) && (expand->elmceq != elmceqRATIO));
772 sprintf(err_buf, "weight-equil-number-all-lambda (%d) must be a positive integer if lmc-weights-equil=%s",
773 expand->equil_n_at_lam, elmceq_names[elmceqNUMATLAM]);
774 CHECK((expand->equil_n_at_lam <= 0) && (expand->elmceq == elmceqNUMATLAM));
776 sprintf(err_buf, "weight-equil-number-samples (%d) must be a positive integer if lmc-weights-equil=%s",
777 expand->equil_samples, elmceq_names[elmceqSAMPLES]);
778 CHECK((expand->equil_samples <= 0) && (expand->elmceq == elmceqSAMPLES));
780 sprintf(err_buf, "weight-equil-number-steps (%d) must be a positive integer if lmc-weights-equil=%s",
781 expand->equil_steps, elmceq_names[elmceqSTEPS]);
782 CHECK((expand->equil_steps <= 0) && (expand->elmceq == elmceqSTEPS));
784 sprintf(err_buf, "weight-equil-wl-delta (%f) must be > 0 if lmc-weights-equil=%s",
785 expand->equil_wl_delta, elmceq_names[elmceqWLDELTA]);
786 CHECK((expand->equil_wl_delta <= 0) && (expand->elmceq == elmceqWLDELTA));
788 sprintf(err_buf, "weight-equil-count-ratio (%f) must be > 0 if lmc-weights-equil=%s",
789 expand->equil_ratio, elmceq_names[elmceqRATIO]);
790 CHECK((expand->equil_ratio <= 0) && (expand->elmceq == elmceqRATIO));
792 sprintf(err_buf, "lmc-weights-equil=%s only possible when lmc-stats = %s or lmc-stats %s",
793 elmceq_names[elmceqWLDELTA], elamstats_names[elamstatsWL], elamstats_names[elamstatsWWL]);
794 CHECK((expand->elmceq == elmceqWLDELTA) && (!EWL(expand->elamstats)));
796 sprintf(err_buf, "lmc-repeats (%d) must be greater than 0", expand->lmc_repeats);
797 CHECK((expand->lmc_repeats <= 0));
798 sprintf(err_buf, "minimum-var-min (%d) must be greater than 0", expand->minvarmin);
799 CHECK((expand->minvarmin <= 0));
800 sprintf(err_buf, "weight-c-range (%d) must be greater or equal to 0", expand->c_range);
801 CHECK((expand->c_range < 0));
802 sprintf(err_buf, "init-lambda-state (%d) must be zero if lmc-forced-nstart (%d)> 0 and lmc-move != 'no'",
803 fep->init_fep_state, expand->lmc_forced_nstart);
804 CHECK((fep->init_fep_state != 0) && (expand->lmc_forced_nstart > 0) && (expand->elmcmove != elmcmoveNO));
805 sprintf(err_buf, "lmc-forced-nstart (%d) must not be negative", expand->lmc_forced_nstart);
806 CHECK((expand->lmc_forced_nstart < 0));
807 sprintf(err_buf, "init-lambda-state (%d) must be in the interval [0,number of lambdas)", fep->init_fep_state);
808 CHECK((fep->init_fep_state < 0) || (fep->init_fep_state >= fep->n_lambda));
810 sprintf(err_buf, "init-wl-delta (%f) must be greater than or equal to 0", expand->init_wl_delta);
811 CHECK((expand->init_wl_delta < 0));
812 sprintf(err_buf, "wl-ratio (%f) must be between 0 and 1", expand->wl_ratio);
813 CHECK((expand->wl_ratio <= 0) || (expand->wl_ratio >= 1));
814 sprintf(err_buf, "wl-scale (%f) must be between 0 and 1", expand->wl_scale);
815 CHECK((expand->wl_scale <= 0) || (expand->wl_scale >= 1));
817 /* if there is no temperature control, we need to specify an MC temperature */
818 sprintf(err_buf, "If there is no temperature control, and lmc-mcmove!= 'no',mc_temperature must be set to a positive number");
819 if (expand->nstTij > 0)
821 sprintf(err_buf, "nst-transition-matrix (%d) must be an integer multiple of nstlog (%d)",
822 expand->nstTij, ir->nstlog);
823 CHECK((mod(expand->nstTij, ir->nstlog) != 0));
828 sprintf(err_buf, "walls only work with pbc=%s", epbc_names[epbcXY]);
829 CHECK(ir->nwall && ir->ePBC != epbcXY);
832 if (ir->ePBC != epbcXYZ && ir->nwall != 2)
834 if (ir->ePBC == epbcNONE)
836 if (ir->epc != epcNO)
838 warning(wi, "Turning off pressure coupling for vacuum system");
844 sprintf(err_buf, "Can not have pressure coupling with pbc=%s",
845 epbc_names[ir->ePBC]);
846 CHECK(ir->epc != epcNO);
848 sprintf(err_buf, "Can not have Ewald with pbc=%s", epbc_names[ir->ePBC]);
849 CHECK(EEL_FULL(ir->coulombtype));
851 sprintf(err_buf, "Can not have dispersion correction with pbc=%s",
852 epbc_names[ir->ePBC]);
853 CHECK(ir->eDispCorr != edispcNO);
856 if (ir->rlist == 0.0)
858 sprintf(err_buf, "can only have neighborlist cut-off zero (=infinite)\n"
859 "with coulombtype = %s or coulombtype = %s\n"
860 "without periodic boundary conditions (pbc = %s) and\n"
861 "rcoulomb and rvdw set to zero",
862 eel_names[eelCUT], eel_names[eelUSER], epbc_names[epbcNONE]);
863 CHECK(((ir->coulombtype != eelCUT) && (ir->coulombtype != eelUSER)) ||
864 (ir->ePBC != epbcNONE) ||
865 (ir->rcoulomb != 0.0) || (ir->rvdw != 0.0));
869 warning_error(wi, "Can not have heuristic neighborlist updates without cut-off");
873 warning_note(wi, "Simulating without cut-offs can be (slightly) faster with nstlist=0, nstype=simple and only one MPI rank");
878 if (ir->nstcomm == 0)
880 ir->comm_mode = ecmNO;
882 if (ir->comm_mode != ecmNO)
886 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");
887 ir->nstcomm = abs(ir->nstcomm);
890 if (ir->nstcalcenergy > 0 && ir->nstcomm < ir->nstcalcenergy)
892 warning_note(wi, "nstcomm < nstcalcenergy defeats the purpose of nstcalcenergy, setting nstcomm to nstcalcenergy");
893 ir->nstcomm = ir->nstcalcenergy;
896 if (ir->comm_mode == ecmANGULAR)
898 sprintf(err_buf, "Can not remove the rotation around the center of mass with periodic molecules");
899 CHECK(ir->bPeriodicMols);
900 if (ir->ePBC != epbcNONE)
902 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).");
907 if (EI_STATE_VELOCITY(ir->eI) && ir->ePBC == epbcNONE && ir->comm_mode != ecmANGULAR)
909 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.");
912 sprintf(err_buf, "Twin-range neighbour searching (NS) with simple NS"
913 " algorithm not implemented");
914 CHECK(((ir->rcoulomb > ir->rlist) || (ir->rvdw > ir->rlist))
915 && (ir->ns_type == ensSIMPLE));
917 /* TEMPERATURE COUPLING */
918 if (ir->etc == etcYES)
920 ir->etc = etcBERENDSEN;
921 warning_note(wi, "Old option for temperature coupling given: "
922 "changing \"yes\" to \"Berendsen\"\n");
925 if ((ir->etc == etcNOSEHOOVER) || (ir->epc == epcMTTK))
927 if (ir->opts.nhchainlength < 1)
929 sprintf(warn_buf, "number of Nose-Hoover chains (currently %d) cannot be less than 1,reset to 1\n", ir->opts.nhchainlength);
930 ir->opts.nhchainlength = 1;
931 warning(wi, warn_buf);
934 if (ir->etc == etcNOSEHOOVER && !EI_VV(ir->eI) && ir->opts.nhchainlength > 1)
936 warning_note(wi, "leapfrog does not yet support Nose-Hoover chains, nhchainlength reset to 1");
937 ir->opts.nhchainlength = 1;
942 ir->opts.nhchainlength = 0;
945 if (ir->eI == eiVVAK)
947 sprintf(err_buf, "%s implemented primarily for validation, and requires nsttcouple = 1 and nstpcouple = 1.",
949 CHECK((ir->nsttcouple != 1) || (ir->nstpcouple != 1));
952 if (ETC_ANDERSEN(ir->etc))
954 sprintf(err_buf, "%s temperature control not supported for integrator %s.", etcoupl_names[ir->etc], ei_names[ir->eI]);
955 CHECK(!(EI_VV(ir->eI)));
957 for (i = 0; i < ir->opts.ngtc; i++)
959 sprintf(err_buf, "all tau_t must currently be equal using Andersen temperature control, violated for group %d", i);
960 CHECK(ir->opts.tau_t[0] != ir->opts.tau_t[i]);
961 sprintf(err_buf, "all tau_t must be postive using Andersen temperature control, tau_t[%d]=%10.6f",
962 i, ir->opts.tau_t[i]);
963 CHECK(ir->opts.tau_t[i] < 0);
965 if (ir->nstcomm > 0 && (ir->etc == etcANDERSEN))
967 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]);
968 warning_note(wi, warn_buf);
971 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]);
972 CHECK(ir->nstcomm > 1 && (ir->etc == etcANDERSEN));
974 for (i = 0; i < ir->opts.ngtc; i++)
976 int nsteps = (int)(ir->opts.tau_t[i]/ir->delta_t);
977 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);
978 CHECK((nsteps % ir->nstcomm) && (ir->etc == etcANDERSENMASSIVE));
981 if (ir->etc == etcBERENDSEN)
983 sprintf(warn_buf, "The %s thermostat does not generate the correct kinetic energy distribution. You might want to consider using the %s thermostat.",
984 ETCOUPLTYPE(ir->etc), ETCOUPLTYPE(etcVRESCALE));
985 warning_note(wi, warn_buf);
988 if ((ir->etc == etcNOSEHOOVER || ETC_ANDERSEN(ir->etc))
989 && ir->epc == epcBERENDSEN)
991 sprintf(warn_buf, "Using Berendsen pressure coupling invalidates the "
992 "true ensemble for the thermostat");
993 warning(wi, warn_buf);
996 /* PRESSURE COUPLING */
997 if (ir->epc == epcISOTROPIC)
999 ir->epc = epcBERENDSEN;
1000 warning_note(wi, "Old option for pressure coupling given: "
1001 "changing \"Isotropic\" to \"Berendsen\"\n");
1004 if (ir->epc != epcNO)
1006 dt_pcoupl = ir->nstpcouple*ir->delta_t;
1008 sprintf(err_buf, "tau-p must be > 0 instead of %g\n", ir->tau_p);
1009 CHECK(ir->tau_p <= 0);
1011 if (ir->tau_p/dt_pcoupl < pcouple_min_integration_steps(ir->epc))
1013 sprintf(warn_buf, "For proper integration of the %s barostat, tau-p (%g) should be at least %d times larger than nstpcouple*dt (%g)",
1014 EPCOUPLTYPE(ir->epc), ir->tau_p, pcouple_min_integration_steps(ir->epc), dt_pcoupl);
1015 warning(wi, warn_buf);
1018 sprintf(err_buf, "compressibility must be > 0 when using pressure"
1019 " coupling %s\n", EPCOUPLTYPE(ir->epc));
1020 CHECK(ir->compress[XX][XX] < 0 || ir->compress[YY][YY] < 0 ||
1021 ir->compress[ZZ][ZZ] < 0 ||
1022 (trace(ir->compress) == 0 && ir->compress[YY][XX] <= 0 &&
1023 ir->compress[ZZ][XX] <= 0 && ir->compress[ZZ][YY] <= 0));
1025 if (epcPARRINELLORAHMAN == ir->epc && opts->bGenVel)
1028 "You are generating velocities so I am assuming you "
1029 "are equilibrating a system. You are using "
1030 "%s pressure coupling, but this can be "
1031 "unstable for equilibration. If your system crashes, try "
1032 "equilibrating first with Berendsen pressure coupling. If "
1033 "you are not equilibrating the system, you can probably "
1034 "ignore this warning.",
1035 epcoupl_names[ir->epc]);
1036 warning(wi, warn_buf);
1042 if (ir->epc > epcNO)
1044 if ((ir->epc != epcBERENDSEN) && (ir->epc != epcMTTK))
1046 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.");
1051 /* ELECTROSTATICS */
1052 /* More checks are in triple check (grompp.c) */
1054 if (ir->coulombtype == eelSWITCH)
1056 sprintf(warn_buf, "coulombtype = %s is only for testing purposes and can lead to serious "
1057 "artifacts, advice: use coulombtype = %s",
1058 eel_names[ir->coulombtype],
1059 eel_names[eelRF_ZERO]);
1060 warning(wi, warn_buf);
1063 if (ir->epsilon_r != 1 && ir->implicit_solvent == eisGBSA)
1065 sprintf(warn_buf, "epsilon-r = %g with GB implicit solvent, will use this value for inner dielectric", ir->epsilon_r);
1066 warning_note(wi, warn_buf);
1069 if (EEL_RF(ir->coulombtype) && ir->epsilon_rf == 1 && ir->epsilon_r != 1)
1071 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);
1072 warning(wi, warn_buf);
1073 ir->epsilon_rf = ir->epsilon_r;
1074 ir->epsilon_r = 1.0;
1077 if (getenv("GALACTIC_DYNAMICS") == NULL)
1079 sprintf(err_buf, "epsilon-r must be >= 0 instead of %g\n", ir->epsilon_r);
1080 CHECK(ir->epsilon_r < 0);
1083 if (EEL_RF(ir->coulombtype))
1085 /* reaction field (at the cut-off) */
1087 if (ir->coulombtype == eelRF_ZERO)
1089 sprintf(warn_buf, "With coulombtype = %s, epsilon-rf must be 0, assuming you meant epsilon_rf=0",
1090 eel_names[ir->coulombtype]);
1091 CHECK(ir->epsilon_rf != 0);
1092 ir->epsilon_rf = 0.0;
1095 sprintf(err_buf, "epsilon-rf must be >= epsilon-r");
1096 CHECK((ir->epsilon_rf < ir->epsilon_r && ir->epsilon_rf != 0) ||
1097 (ir->epsilon_r == 0));
1098 if (ir->epsilon_rf == ir->epsilon_r)
1100 sprintf(warn_buf, "Using epsilon-rf = epsilon-r with %s does not make sense",
1101 eel_names[ir->coulombtype]);
1102 warning(wi, warn_buf);
1105 /* Allow rlist>rcoulomb for tabulated long range stuff. This just
1106 * means the interaction is zero outside rcoulomb, but it helps to
1107 * provide accurate energy conservation.
1109 if (ir_coulomb_might_be_zero_at_cutoff(ir))
1111 if (ir_coulomb_switched(ir))
1114 "With coulombtype = %s rcoulomb_switch must be < rcoulomb. Or, better: Use the potential modifier options!",
1115 eel_names[ir->coulombtype]);
1116 CHECK(ir->rcoulomb_switch >= ir->rcoulomb);
1119 else if (ir->coulombtype == eelCUT || EEL_RF(ir->coulombtype))
1121 if (ir->cutoff_scheme == ecutsGROUP && ir->coulomb_modifier == eintmodNONE)
1123 sprintf(err_buf, "With coulombtype = %s, rcoulomb should be >= rlist unless you use a potential modifier",
1124 eel_names[ir->coulombtype]);
1125 CHECK(ir->rlist > ir->rcoulomb);
1129 if (ir->coulombtype == eelSWITCH || ir->coulombtype == eelSHIFT ||
1130 ir->vdwtype == evdwSWITCH || ir->vdwtype == evdwSHIFT)
1133 "The switch/shift interaction settings are just for compatibility; you will get better "
1134 "performance from applying potential modifiers to your interactions!\n");
1135 warning_note(wi, warn_buf);
1138 if (ir->coulombtype == eelPMESWITCH)
1140 if (ir->rcoulomb_switch/ir->rcoulomb < 0.9499)
1142 sprintf(warn_buf, "The switching range for %s should be 5%% or less, energy conservation will be good anyhow, since ewald_rtol = %g",
1143 eel_names[ir->coulombtype],
1145 warning(wi, warn_buf);
1149 if (EEL_FULL(ir->coulombtype))
1151 if (ir->coulombtype == eelPMESWITCH || ir->coulombtype == eelPMEUSER ||
1152 ir->coulombtype == eelPMEUSERSWITCH)
1154 sprintf(err_buf, "With coulombtype = %s, rcoulomb must be <= rlist",
1155 eel_names[ir->coulombtype]);
1156 CHECK(ir->rcoulomb > ir->rlist);
1158 else if (ir->cutoff_scheme == ecutsGROUP && ir->coulomb_modifier == eintmodNONE)
1160 if (ir->coulombtype == eelPME || ir->coulombtype == eelP3M_AD)
1163 "With coulombtype = %s (without modifier), rcoulomb must be equal to rlist,\n"
1164 "or rlistlong if nstcalclr=1. For optimal energy conservation,consider using\n"
1165 "a potential modifier.", eel_names[ir->coulombtype]);
1166 if (ir->nstcalclr == 1)
1168 CHECK(ir->rcoulomb != ir->rlist && ir->rcoulomb != ir->rlistlong);
1172 CHECK(ir->rcoulomb != ir->rlist);
1178 if (EEL_PME(ir->coulombtype) || EVDW_PME(ir->vdwtype))
1180 if (ir->pme_order < 3)
1182 warning_error(wi, "pme-order can not be smaller than 3");
1186 if (ir->nwall == 2 && EEL_FULL(ir->coulombtype))
1188 if (ir->ewald_geometry == eewg3D)
1190 sprintf(warn_buf, "With pbc=%s you should use ewald-geometry=%s",
1191 epbc_names[ir->ePBC], eewg_names[eewg3DC]);
1192 warning(wi, warn_buf);
1194 /* This check avoids extra pbc coding for exclusion corrections */
1195 sprintf(err_buf, "wall-ewald-zfac should be >= 2");
1196 CHECK(ir->wall_ewald_zfac < 2);
1199 if (ir_vdw_switched(ir))
1201 sprintf(err_buf, "With switched vdw forces or potentials, rvdw-switch must be < rvdw");
1202 CHECK(ir->rvdw_switch >= ir->rvdw);
1204 if (ir->rvdw_switch < 0.5*ir->rvdw)
1206 sprintf(warn_buf, "You are applying a switch function to vdw forces or potentials from %g to %g nm, which is more than half the interaction range, whereas switch functions are intended to act only close to the cut-off.",
1207 ir->rvdw_switch, ir->rvdw);
1208 warning_note(wi, warn_buf);
1211 else if (ir->vdwtype == evdwCUT || ir->vdwtype == evdwPME)
1213 if (ir->cutoff_scheme == ecutsGROUP && ir->vdw_modifier == eintmodNONE)
1215 sprintf(err_buf, "With vdwtype = %s, rvdw must be >= rlist unless you use a potential modifier", evdw_names[ir->vdwtype]);
1216 CHECK(ir->rlist > ir->rvdw);
1220 if (ir->cutoff_scheme == ecutsGROUP)
1222 if (((ir->coulomb_modifier != eintmodNONE && ir->rcoulomb == ir->rlist) ||
1223 (ir->vdw_modifier != eintmodNONE && ir->rvdw == ir->rlist)) &&
1226 warning_note(wi, "With exact cut-offs, rlist should be "
1227 "larger than rcoulomb and rvdw, so that there "
1228 "is a buffer region for particle motion "
1229 "between neighborsearch steps");
1232 if (ir_coulomb_is_zero_at_cutoff(ir) && ir->rlistlong <= ir->rcoulomb)
1234 sprintf(warn_buf, "For energy conservation with switch/shift potentials, %s should be 0.1 to 0.3 nm larger than rcoulomb.",
1235 IR_TWINRANGE(*ir) ? "rlistlong" : "rlist");
1236 warning_note(wi, warn_buf);
1238 if (ir_vdw_switched(ir) && (ir->rlistlong <= ir->rvdw))
1240 sprintf(warn_buf, "For energy conservation with switch/shift potentials, %s should be 0.1 to 0.3 nm larger than rvdw.",
1241 IR_TWINRANGE(*ir) ? "rlistlong" : "rlist");
1242 warning_note(wi, warn_buf);
1246 if (ir->vdwtype == evdwUSER && ir->eDispCorr != edispcNO)
1248 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.");
1251 if (ir->nstlist == -1)
1253 sprintf(err_buf, "With nstlist=-1 rvdw and rcoulomb should be smaller than rlist to account for diffusion and possibly charge-group radii");
1254 CHECK(ir->rvdw >= ir->rlist || ir->rcoulomb >= ir->rlist);
1256 sprintf(err_buf, "nstlist can not be smaller than -1");
1257 CHECK(ir->nstlist < -1);
1259 if (ir->eI == eiLBFGS && (ir->coulombtype == eelCUT || ir->vdwtype == evdwCUT)
1262 warning(wi, "For efficient BFGS minimization, use switch/shift/pme instead of cut-off.");
1265 if (ir->eI == eiLBFGS && ir->nbfgscorr <= 0)
1267 warning(wi, "Using L-BFGS with nbfgscorr<=0 just gets you steepest descent.");
1270 /* ENERGY CONSERVATION */
1271 if (ir_NVE(ir) && ir->cutoff_scheme == ecutsGROUP)
1273 if (!ir_vdw_might_be_zero_at_cutoff(ir) && ir->rvdw > 0 && ir->vdw_modifier == eintmodNONE)
1275 sprintf(warn_buf, "You are using a cut-off for VdW interactions with NVE, for good energy conservation use vdwtype = %s (possibly with DispCorr)",
1276 evdw_names[evdwSHIFT]);
1277 warning_note(wi, warn_buf);
1279 if (!ir_coulomb_might_be_zero_at_cutoff(ir) && ir->rcoulomb > 0)
1281 sprintf(warn_buf, "You are using a cut-off for electrostatics with NVE, for good energy conservation use coulombtype = %s or %s",
1282 eel_names[eelPMESWITCH], eel_names[eelRF_ZERO]);
1283 warning_note(wi, warn_buf);
1287 /* IMPLICIT SOLVENT */
1288 if (ir->coulombtype == eelGB_NOTUSED)
1290 ir->coulombtype = eelCUT;
1291 ir->implicit_solvent = eisGBSA;
1292 fprintf(stderr, "Note: Old option for generalized born electrostatics given:\n"
1293 "Changing coulombtype from \"generalized-born\" to \"cut-off\" and instead\n"
1294 "setting implicit-solvent value to \"GBSA\" in input section.\n");
1297 if (ir->sa_algorithm == esaSTILL)
1299 sprintf(err_buf, "Still SA algorithm not available yet, use %s or %s instead\n", esa_names[esaAPPROX], esa_names[esaNO]);
1300 CHECK(ir->sa_algorithm == esaSTILL);
1303 if (ir->implicit_solvent == eisGBSA)
1305 sprintf(err_buf, "With GBSA implicit solvent, rgbradii must be equal to rlist.");
1306 CHECK(ir->rgbradii != ir->rlist);
1308 if (ir->coulombtype != eelCUT)
1310 sprintf(err_buf, "With GBSA, coulombtype must be equal to %s\n", eel_names[eelCUT]);
1311 CHECK(ir->coulombtype != eelCUT);
1313 if (ir->vdwtype != evdwCUT)
1315 sprintf(err_buf, "With GBSA, vdw-type must be equal to %s\n", evdw_names[evdwCUT]);
1316 CHECK(ir->vdwtype != evdwCUT);
1318 if (ir->nstgbradii < 1)
1320 sprintf(warn_buf, "Using GBSA with nstgbradii<1, setting nstgbradii=1");
1321 warning_note(wi, warn_buf);
1324 if (ir->sa_algorithm == esaNO)
1326 sprintf(warn_buf, "No SA (non-polar) calculation requested together with GB. Are you sure this is what you want?\n");
1327 warning_note(wi, warn_buf);
1329 if (ir->sa_surface_tension < 0 && ir->sa_algorithm != esaNO)
1331 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");
1332 warning_note(wi, warn_buf);
1334 if (ir->gb_algorithm == egbSTILL)
1336 ir->sa_surface_tension = 0.0049 * CAL2JOULE * 100;
1340 ir->sa_surface_tension = 0.0054 * CAL2JOULE * 100;
1343 if (ir->sa_surface_tension == 0 && ir->sa_algorithm != esaNO)
1345 sprintf(err_buf, "Surface tension set to 0 while SA-calculation requested\n");
1346 CHECK(ir->sa_surface_tension == 0 && ir->sa_algorithm != esaNO);
1353 if (ir->cutoff_scheme != ecutsGROUP)
1355 warning_error(wi, "AdresS simulation supports only cutoff-scheme=group");
1359 warning_error(wi, "AdresS simulation supports only stochastic dynamics");
1361 if (ir->epc != epcNO)
1363 warning_error(wi, "AdresS simulation does not support pressure coupling");
1365 if (EEL_FULL(ir->coulombtype))
1367 warning_error(wi, "AdresS simulation does not support long-range electrostatics");
1372 /* count the number of text elemets separated by whitespace in a string.
1373 str = the input string
1374 maxptr = the maximum number of allowed elements
1375 ptr = the output array of pointers to the first character of each element
1376 returns: the number of elements. */
1377 int str_nelem(const char *str, int maxptr, char *ptr[])
1382 copy0 = strdup(str);
1385 while (*copy != '\0')
1389 gmx_fatal(FARGS, "Too many groups on line: '%s' (max is %d)",
1397 while ((*copy != '\0') && !isspace(*copy))
1416 /* interpret a number of doubles from a string and put them in an array,
1417 after allocating space for them.
1418 str = the input string
1419 n = the (pre-allocated) number of doubles read
1420 r = the output array of doubles. */
1421 static void parse_n_real(char *str, int *n, real **r)
1426 *n = str_nelem(str, MAXPTR, ptr);
1429 for (i = 0; i < *n; i++)
1431 (*r)[i] = strtod(ptr[i], NULL);
1435 static void do_fep_params(t_inputrec *ir, char fep_lambda[][STRLEN], char weights[STRLEN])
1438 int i, j, max_n_lambda, nweights, nfep[efptNR];
1439 t_lambda *fep = ir->fepvals;
1440 t_expanded *expand = ir->expandedvals;
1441 real **count_fep_lambdas;
1442 gmx_bool bOneLambda = TRUE;
1444 snew(count_fep_lambdas, efptNR);
1446 /* FEP input processing */
1447 /* first, identify the number of lambda values for each type.
1448 All that are nonzero must have the same number */
1450 for (i = 0; i < efptNR; i++)
1452 parse_n_real(fep_lambda[i], &(nfep[i]), &(count_fep_lambdas[i]));
1455 /* now, determine the number of components. All must be either zero, or equal. */
1458 for (i = 0; i < efptNR; i++)
1460 if (nfep[i] > max_n_lambda)
1462 max_n_lambda = nfep[i]; /* here's a nonzero one. All of them
1463 must have the same number if its not zero.*/
1468 for (i = 0; i < efptNR; i++)
1472 ir->fepvals->separate_dvdl[i] = FALSE;
1474 else if (nfep[i] == max_n_lambda)
1476 if (i != efptTEMPERATURE) /* we treat this differently -- not really a reason to compute the derivative with
1477 respect to the temperature currently */
1479 ir->fepvals->separate_dvdl[i] = TRUE;
1484 gmx_fatal(FARGS, "Number of lambdas (%d) for FEP type %s not equal to number of other types (%d)",
1485 nfep[i], efpt_names[i], max_n_lambda);
1488 /* we don't print out dhdl if the temperature is changing, since we can't correctly define dhdl in this case */
1489 ir->fepvals->separate_dvdl[efptTEMPERATURE] = FALSE;
1491 /* the number of lambdas is the number we've read in, which is either zero
1492 or the same for all */
1493 fep->n_lambda = max_n_lambda;
1495 /* allocate space for the array of lambda values */
1496 snew(fep->all_lambda, efptNR);
1497 /* if init_lambda is defined, we need to set lambda */
1498 if ((fep->init_lambda > 0) && (fep->n_lambda == 0))
1500 ir->fepvals->separate_dvdl[efptFEP] = TRUE;
1502 /* otherwise allocate the space for all of the lambdas, and transfer the data */
1503 for (i = 0; i < efptNR; i++)
1505 snew(fep->all_lambda[i], fep->n_lambda);
1506 if (nfep[i] > 0) /* if it's zero, then the count_fep_lambda arrays
1509 for (j = 0; j < fep->n_lambda; j++)
1511 fep->all_lambda[i][j] = (double)count_fep_lambdas[i][j];
1513 sfree(count_fep_lambdas[i]);
1516 sfree(count_fep_lambdas);
1518 /* "fep-vals" is either zero or the full number. If zero, we'll need to define fep-lambdas for internal
1519 bookkeeping -- for now, init_lambda */
1521 if ((nfep[efptFEP] == 0) && (fep->init_lambda >= 0))
1523 for (i = 0; i < fep->n_lambda; i++)
1525 fep->all_lambda[efptFEP][i] = fep->init_lambda;
1529 /* check to see if only a single component lambda is defined, and soft core is defined.
1530 In this case, turn on coulomb soft core */
1532 if (max_n_lambda == 0)
1538 for (i = 0; i < efptNR; i++)
1540 if ((nfep[i] != 0) && (i != efptFEP))
1546 if ((bOneLambda) && (fep->sc_alpha > 0))
1548 fep->bScCoul = TRUE;
1551 /* Fill in the others with the efptFEP if they are not explicitly
1552 specified (i.e. nfep[i] == 0). This means if fep is not defined,
1553 they are all zero. */
1555 for (i = 0; i < efptNR; i++)
1557 if ((nfep[i] == 0) && (i != efptFEP))
1559 for (j = 0; j < fep->n_lambda; j++)
1561 fep->all_lambda[i][j] = fep->all_lambda[efptFEP][j];
1567 /* make it easier if sc_r_power = 48 by increasing it to the 4th power, to be in the right scale. */
1568 if (fep->sc_r_power == 48)
1570 if (fep->sc_alpha > 0.1)
1572 gmx_fatal(FARGS, "sc_alpha (%f) for sc_r_power = 48 should usually be between 0.001 and 0.004", fep->sc_alpha);
1576 expand = ir->expandedvals;
1577 /* now read in the weights */
1578 parse_n_real(weights, &nweights, &(expand->init_lambda_weights));
1581 snew(expand->init_lambda_weights, fep->n_lambda); /* initialize to zero */
1583 else if (nweights != fep->n_lambda)
1585 gmx_fatal(FARGS, "Number of weights (%d) is not equal to number of lambda values (%d)",
1586 nweights, fep->n_lambda);
1588 if ((expand->nstexpanded < 0) && (ir->efep != efepNO))
1590 expand->nstexpanded = fep->nstdhdl;
1591 /* if you don't specify nstexpanded when doing expanded ensemble free energy calcs, it is set to nstdhdl */
1593 if ((expand->nstexpanded < 0) && ir->bSimTemp)
1595 expand->nstexpanded = 2*(int)(ir->opts.tau_t[0]/ir->delta_t);
1596 /* if you don't specify nstexpanded when doing expanded ensemble simulated tempering, it is set to
1597 2*tau_t just to be careful so it's not to frequent */
1602 static void do_simtemp_params(t_inputrec *ir)
1605 snew(ir->simtempvals->temperatures, ir->fepvals->n_lambda);
1606 GetSimTemps(ir->fepvals->n_lambda, ir->simtempvals, ir->fepvals->all_lambda[efptTEMPERATURE]);
1611 static void do_wall_params(t_inputrec *ir,
1612 char *wall_atomtype, char *wall_density,
1616 char *names[MAXPTR];
1619 opts->wall_atomtype[0] = NULL;
1620 opts->wall_atomtype[1] = NULL;
1622 ir->wall_atomtype[0] = -1;
1623 ir->wall_atomtype[1] = -1;
1624 ir->wall_density[0] = 0;
1625 ir->wall_density[1] = 0;
1629 nstr = str_nelem(wall_atomtype, MAXPTR, names);
1630 if (nstr != ir->nwall)
1632 gmx_fatal(FARGS, "Expected %d elements for wall_atomtype, found %d",
1635 for (i = 0; i < ir->nwall; i++)
1637 opts->wall_atomtype[i] = strdup(names[i]);
1640 if (ir->wall_type == ewt93 || ir->wall_type == ewt104)
1642 nstr = str_nelem(wall_density, MAXPTR, names);
1643 if (nstr != ir->nwall)
1645 gmx_fatal(FARGS, "Expected %d elements for wall-density, found %d", ir->nwall, nstr);
1647 for (i = 0; i < ir->nwall; i++)
1649 sscanf(names[i], "%lf", &dbl);
1652 gmx_fatal(FARGS, "wall-density[%d] = %f\n", i, dbl);
1654 ir->wall_density[i] = dbl;
1660 static void add_wall_energrps(gmx_groups_t *groups, int nwall, t_symtab *symtab)
1668 srenew(groups->grpname, groups->ngrpname+nwall);
1669 grps = &(groups->grps[egcENER]);
1670 srenew(grps->nm_ind, grps->nr+nwall);
1671 for (i = 0; i < nwall; i++)
1673 sprintf(str, "wall%d", i);
1674 groups->grpname[groups->ngrpname] = put_symtab(symtab, str);
1675 grps->nm_ind[grps->nr++] = groups->ngrpname++;
1680 void read_expandedparams(int *ninp_p, t_inpfile **inp_p,
1681 t_expanded *expand, warninp_t wi)
1683 int ninp, nerror = 0;
1689 /* read expanded ensemble parameters */
1690 CCTYPE ("expanded ensemble variables");
1691 ITYPE ("nstexpanded", expand->nstexpanded, -1);
1692 EETYPE("lmc-stats", expand->elamstats, elamstats_names);
1693 EETYPE("lmc-move", expand->elmcmove, elmcmove_names);
1694 EETYPE("lmc-weights-equil", expand->elmceq, elmceq_names);
1695 ITYPE ("weight-equil-number-all-lambda", expand->equil_n_at_lam, -1);
1696 ITYPE ("weight-equil-number-samples", expand->equil_samples, -1);
1697 ITYPE ("weight-equil-number-steps", expand->equil_steps, -1);
1698 RTYPE ("weight-equil-wl-delta", expand->equil_wl_delta, -1);
1699 RTYPE ("weight-equil-count-ratio", expand->equil_ratio, -1);
1700 CCTYPE("Seed for Monte Carlo in lambda space");
1701 ITYPE ("lmc-seed", expand->lmc_seed, -1);
1702 RTYPE ("mc-temperature", expand->mc_temp, -1);
1703 ITYPE ("lmc-repeats", expand->lmc_repeats, 1);
1704 ITYPE ("lmc-gibbsdelta", expand->gibbsdeltalam, -1);
1705 ITYPE ("lmc-forced-nstart", expand->lmc_forced_nstart, 0);
1706 EETYPE("symmetrized-transition-matrix", expand->bSymmetrizedTMatrix, yesno_names);
1707 ITYPE("nst-transition-matrix", expand->nstTij, -1);
1708 ITYPE ("mininum-var-min", expand->minvarmin, 100); /*default is reasonable */
1709 ITYPE ("weight-c-range", expand->c_range, 0); /* default is just C=0 */
1710 RTYPE ("wl-scale", expand->wl_scale, 0.8);
1711 RTYPE ("wl-ratio", expand->wl_ratio, 0.8);
1712 RTYPE ("init-wl-delta", expand->init_wl_delta, 1.0);
1713 EETYPE("wl-oneovert", expand->bWLoneovert, yesno_names);
1721 void get_ir(const char *mdparin, const char *mdparout,
1722 t_inputrec *ir, t_gromppopts *opts,
1726 double dumdub[2][6];
1730 char warn_buf[STRLEN];
1731 t_lambda *fep = ir->fepvals;
1732 t_expanded *expand = ir->expandedvals;
1734 init_inputrec_strings();
1735 inp = read_inpfile(mdparin, &ninp, wi);
1737 snew(dumstr[0], STRLEN);
1738 snew(dumstr[1], STRLEN);
1740 if (-1 == search_einp(ninp, inp, "cutoff-scheme"))
1743 "%s did not specify a value for the .mdp option "
1744 "\"cutoff-scheme\". Probably it was first intended for use "
1745 "with GROMACS before 4.6. In 4.6, the Verlet scheme was "
1746 "introduced, but the group scheme was still the default. "
1747 "The default is now the Verlet scheme, so you will observe "
1748 "different behaviour.", mdparin);
1749 warning_note(wi, warn_buf);
1752 /* remove the following deprecated commands */
1755 REM_TYPE("domain-decomposition");
1756 REM_TYPE("andersen-seed");
1758 REM_TYPE("dihre-fc");
1759 REM_TYPE("dihre-tau");
1760 REM_TYPE("nstdihreout");
1761 REM_TYPE("nstcheckpoint");
1763 /* replace the following commands with the clearer new versions*/
1764 REPL_TYPE("unconstrained-start", "continuation");
1765 REPL_TYPE("foreign-lambda", "fep-lambdas");
1766 REPL_TYPE("verlet-buffer-drift", "verlet-buffer-tolerance");
1767 REPL_TYPE("nstxtcout", "nstxout-compressed");
1768 REPL_TYPE("xtc-grps", "compressed-x-grps");
1769 REPL_TYPE("xtc-precision", "compressed-x-precision");
1771 CCTYPE ("VARIOUS PREPROCESSING OPTIONS");
1772 CTYPE ("Preprocessor information: use cpp syntax.");
1773 CTYPE ("e.g.: -I/home/joe/doe -I/home/mary/roe");
1774 STYPE ("include", opts->include, NULL);
1775 CTYPE ("e.g.: -DPOSRES -DFLEXIBLE (note these variable names are case sensitive)");
1776 STYPE ("define", opts->define, NULL);
1778 CCTYPE ("RUN CONTROL PARAMETERS");
1779 EETYPE("integrator", ir->eI, ei_names);
1780 CTYPE ("Start time and timestep in ps");
1781 RTYPE ("tinit", ir->init_t, 0.0);
1782 RTYPE ("dt", ir->delta_t, 0.001);
1783 STEPTYPE ("nsteps", ir->nsteps, 0);
1784 CTYPE ("For exact run continuation or redoing part of a run");
1785 STEPTYPE ("init-step", ir->init_step, 0);
1786 CTYPE ("Part index is updated automatically on checkpointing (keeps files separate)");
1787 ITYPE ("simulation-part", ir->simulation_part, 1);
1788 CTYPE ("mode for center of mass motion removal");
1789 EETYPE("comm-mode", ir->comm_mode, ecm_names);
1790 CTYPE ("number of steps for center of mass motion removal");
1791 ITYPE ("nstcomm", ir->nstcomm, 100);
1792 CTYPE ("group(s) for center of mass motion removal");
1793 STYPE ("comm-grps", is->vcm, NULL);
1795 CCTYPE ("LANGEVIN DYNAMICS OPTIONS");
1796 CTYPE ("Friction coefficient (amu/ps) and random seed");
1797 RTYPE ("bd-fric", ir->bd_fric, 0.0);
1798 STEPTYPE ("ld-seed", ir->ld_seed, -1);
1801 CCTYPE ("ENERGY MINIMIZATION OPTIONS");
1802 CTYPE ("Force tolerance and initial step-size");
1803 RTYPE ("emtol", ir->em_tol, 10.0);
1804 RTYPE ("emstep", ir->em_stepsize, 0.01);
1805 CTYPE ("Max number of iterations in relax-shells");
1806 ITYPE ("niter", ir->niter, 20);
1807 CTYPE ("Step size (ps^2) for minimization of flexible constraints");
1808 RTYPE ("fcstep", ir->fc_stepsize, 0);
1809 CTYPE ("Frequency of steepest descents steps when doing CG");
1810 ITYPE ("nstcgsteep", ir->nstcgsteep, 1000);
1811 ITYPE ("nbfgscorr", ir->nbfgscorr, 10);
1813 CCTYPE ("TEST PARTICLE INSERTION OPTIONS");
1814 RTYPE ("rtpi", ir->rtpi, 0.05);
1816 /* Output options */
1817 CCTYPE ("OUTPUT CONTROL OPTIONS");
1818 CTYPE ("Output frequency for coords (x), velocities (v) and forces (f)");
1819 ITYPE ("nstxout", ir->nstxout, 0);
1820 ITYPE ("nstvout", ir->nstvout, 0);
1821 ITYPE ("nstfout", ir->nstfout, 0);
1822 ir->nstcheckpoint = 1000;
1823 CTYPE ("Output frequency for energies to log file and energy file");
1824 ITYPE ("nstlog", ir->nstlog, 1000);
1825 ITYPE ("nstcalcenergy", ir->nstcalcenergy, 100);
1826 ITYPE ("nstenergy", ir->nstenergy, 1000);
1827 CTYPE ("Output frequency and precision for .xtc file");
1828 ITYPE ("nstxout-compressed", ir->nstxout_compressed, 0);
1829 RTYPE ("compressed-x-precision", ir->x_compression_precision, 1000.0);
1830 CTYPE ("This selects the subset of atoms for the compressed");
1831 CTYPE ("trajectory file. You can select multiple groups. By");
1832 CTYPE ("default, all atoms will be written.");
1833 STYPE ("compressed-x-grps", is->x_compressed_groups, NULL);
1834 CTYPE ("Selection of energy groups");
1835 STYPE ("energygrps", is->energy, NULL);
1837 /* Neighbor searching */
1838 CCTYPE ("NEIGHBORSEARCHING PARAMETERS");
1839 CTYPE ("cut-off scheme (Verlet: particle based cut-offs, group: using charge groups)");
1840 EETYPE("cutoff-scheme", ir->cutoff_scheme, ecutscheme_names);
1841 CTYPE ("nblist update frequency");
1842 ITYPE ("nstlist", ir->nstlist, 10);
1843 CTYPE ("ns algorithm (simple or grid)");
1844 EETYPE("ns-type", ir->ns_type, ens_names);
1845 /* set ndelta to the optimal value of 2 */
1847 CTYPE ("Periodic boundary conditions: xyz, no, xy");
1848 EETYPE("pbc", ir->ePBC, epbc_names);
1849 EETYPE("periodic-molecules", ir->bPeriodicMols, yesno_names);
1850 CTYPE ("Allowed energy error due to the Verlet buffer in kJ/mol/ps per atom,");
1851 CTYPE ("a value of -1 means: use rlist");
1852 RTYPE("verlet-buffer-tolerance", ir->verletbuf_tol, 0.005);
1853 CTYPE ("nblist cut-off");
1854 RTYPE ("rlist", ir->rlist, 1.0);
1855 CTYPE ("long-range cut-off for switched potentials");
1856 RTYPE ("rlistlong", ir->rlistlong, -1);
1857 ITYPE ("nstcalclr", ir->nstcalclr, -1);
1859 /* Electrostatics */
1860 CCTYPE ("OPTIONS FOR ELECTROSTATICS AND VDW");
1861 CTYPE ("Method for doing electrostatics");
1862 EETYPE("coulombtype", ir->coulombtype, eel_names);
1863 EETYPE("coulomb-modifier", ir->coulomb_modifier, eintmod_names);
1864 CTYPE ("cut-off lengths");
1865 RTYPE ("rcoulomb-switch", ir->rcoulomb_switch, 0.0);
1866 RTYPE ("rcoulomb", ir->rcoulomb, 1.0);
1867 CTYPE ("Relative dielectric constant for the medium and the reaction field");
1868 RTYPE ("epsilon-r", ir->epsilon_r, 1.0);
1869 RTYPE ("epsilon-rf", ir->epsilon_rf, 0.0);
1870 CTYPE ("Method for doing Van der Waals");
1871 EETYPE("vdw-type", ir->vdwtype, evdw_names);
1872 EETYPE("vdw-modifier", ir->vdw_modifier, eintmod_names);
1873 CTYPE ("cut-off lengths");
1874 RTYPE ("rvdw-switch", ir->rvdw_switch, 0.0);
1875 RTYPE ("rvdw", ir->rvdw, 1.0);
1876 CTYPE ("Apply long range dispersion corrections for Energy and Pressure");
1877 EETYPE("DispCorr", ir->eDispCorr, edispc_names);
1878 CTYPE ("Extension of the potential lookup tables beyond the cut-off");
1879 RTYPE ("table-extension", ir->tabext, 1.0);
1880 CTYPE ("Separate tables between energy group pairs");
1881 STYPE ("energygrp-table", is->egptable, NULL);
1882 CTYPE ("Spacing for the PME/PPPM FFT grid");
1883 RTYPE ("fourierspacing", ir->fourier_spacing, 0.12);
1884 CTYPE ("FFT grid size, when a value is 0 fourierspacing will be used");
1885 ITYPE ("fourier-nx", ir->nkx, 0);
1886 ITYPE ("fourier-ny", ir->nky, 0);
1887 ITYPE ("fourier-nz", ir->nkz, 0);
1888 CTYPE ("EWALD/PME/PPPM parameters");
1889 ITYPE ("pme-order", ir->pme_order, 4);
1890 RTYPE ("ewald-rtol", ir->ewald_rtol, 0.00001);
1891 RTYPE ("ewald-rtol-lj", ir->ewald_rtol_lj, 0.001);
1892 EETYPE("lj-pme-comb-rule", ir->ljpme_combination_rule, eljpme_names);
1893 EETYPE("ewald-geometry", ir->ewald_geometry, eewg_names);
1894 RTYPE ("epsilon-surface", ir->epsilon_surface, 0.0);
1895 EETYPE("optimize-fft", ir->bOptFFT, yesno_names);
1897 CCTYPE("IMPLICIT SOLVENT ALGORITHM");
1898 EETYPE("implicit-solvent", ir->implicit_solvent, eis_names);
1900 CCTYPE ("GENERALIZED BORN ELECTROSTATICS");
1901 CTYPE ("Algorithm for calculating Born radii");
1902 EETYPE("gb-algorithm", ir->gb_algorithm, egb_names);
1903 CTYPE ("Frequency of calculating the Born radii inside rlist");
1904 ITYPE ("nstgbradii", ir->nstgbradii, 1);
1905 CTYPE ("Cutoff for Born radii calculation; the contribution from atoms");
1906 CTYPE ("between rlist and rgbradii is updated every nstlist steps");
1907 RTYPE ("rgbradii", ir->rgbradii, 1.0);
1908 CTYPE ("Dielectric coefficient of the implicit solvent");
1909 RTYPE ("gb-epsilon-solvent", ir->gb_epsilon_solvent, 80.0);
1910 CTYPE ("Salt concentration in M for Generalized Born models");
1911 RTYPE ("gb-saltconc", ir->gb_saltconc, 0.0);
1912 CTYPE ("Scaling factors used in the OBC GB model. Default values are OBC(II)");
1913 RTYPE ("gb-obc-alpha", ir->gb_obc_alpha, 1.0);
1914 RTYPE ("gb-obc-beta", ir->gb_obc_beta, 0.8);
1915 RTYPE ("gb-obc-gamma", ir->gb_obc_gamma, 4.85);
1916 RTYPE ("gb-dielectric-offset", ir->gb_dielectric_offset, 0.009);
1917 EETYPE("sa-algorithm", ir->sa_algorithm, esa_names);
1918 CTYPE ("Surface tension (kJ/mol/nm^2) for the SA (nonpolar surface) part of GBSA");
1919 CTYPE ("The value -1 will set default value for Still/HCT/OBC GB-models.");
1920 RTYPE ("sa-surface-tension", ir->sa_surface_tension, -1);
1922 /* Coupling stuff */
1923 CCTYPE ("OPTIONS FOR WEAK COUPLING ALGORITHMS");
1924 CTYPE ("Temperature coupling");
1925 EETYPE("tcoupl", ir->etc, etcoupl_names);
1926 ITYPE ("nsttcouple", ir->nsttcouple, -1);
1927 ITYPE("nh-chain-length", ir->opts.nhchainlength, 10);
1928 EETYPE("print-nose-hoover-chain-variables", ir->bPrintNHChains, yesno_names);
1929 CTYPE ("Groups to couple separately");
1930 STYPE ("tc-grps", is->tcgrps, NULL);
1931 CTYPE ("Time constant (ps) and reference temperature (K)");
1932 STYPE ("tau-t", is->tau_t, NULL);
1933 STYPE ("ref-t", is->ref_t, NULL);
1934 CTYPE ("pressure coupling");
1935 EETYPE("pcoupl", ir->epc, epcoupl_names);
1936 EETYPE("pcoupltype", ir->epct, epcoupltype_names);
1937 ITYPE ("nstpcouple", ir->nstpcouple, -1);
1938 CTYPE ("Time constant (ps), compressibility (1/bar) and reference P (bar)");
1939 RTYPE ("tau-p", ir->tau_p, 1.0);
1940 STYPE ("compressibility", dumstr[0], NULL);
1941 STYPE ("ref-p", dumstr[1], NULL);
1942 CTYPE ("Scaling of reference coordinates, No, All or COM");
1943 EETYPE ("refcoord-scaling", ir->refcoord_scaling, erefscaling_names);
1946 CCTYPE ("OPTIONS FOR QMMM calculations");
1947 EETYPE("QMMM", ir->bQMMM, yesno_names);
1948 CTYPE ("Groups treated Quantum Mechanically");
1949 STYPE ("QMMM-grps", is->QMMM, NULL);
1950 CTYPE ("QM method");
1951 STYPE("QMmethod", is->QMmethod, NULL);
1952 CTYPE ("QMMM scheme");
1953 EETYPE("QMMMscheme", ir->QMMMscheme, eQMMMscheme_names);
1954 CTYPE ("QM basisset");
1955 STYPE("QMbasis", is->QMbasis, NULL);
1956 CTYPE ("QM charge");
1957 STYPE ("QMcharge", is->QMcharge, NULL);
1958 CTYPE ("QM multiplicity");
1959 STYPE ("QMmult", is->QMmult, NULL);
1960 CTYPE ("Surface Hopping");
1961 STYPE ("SH", is->bSH, NULL);
1962 CTYPE ("CAS space options");
1963 STYPE ("CASorbitals", is->CASorbitals, NULL);
1964 STYPE ("CASelectrons", is->CASelectrons, NULL);
1965 STYPE ("SAon", is->SAon, NULL);
1966 STYPE ("SAoff", is->SAoff, NULL);
1967 STYPE ("SAsteps", is->SAsteps, NULL);
1968 CTYPE ("Scale factor for MM charges");
1969 RTYPE ("MMChargeScaleFactor", ir->scalefactor, 1.0);
1970 CTYPE ("Optimization of QM subsystem");
1971 STYPE ("bOPT", is->bOPT, NULL);
1972 STYPE ("bTS", is->bTS, NULL);
1974 /* Simulated annealing */
1975 CCTYPE("SIMULATED ANNEALING");
1976 CTYPE ("Type of annealing for each temperature group (no/single/periodic)");
1977 STYPE ("annealing", is->anneal, NULL);
1978 CTYPE ("Number of time points to use for specifying annealing in each group");
1979 STYPE ("annealing-npoints", is->anneal_npoints, NULL);
1980 CTYPE ("List of times at the annealing points for each group");
1981 STYPE ("annealing-time", is->anneal_time, NULL);
1982 CTYPE ("Temp. at each annealing point, for each group.");
1983 STYPE ("annealing-temp", is->anneal_temp, NULL);
1986 CCTYPE ("GENERATE VELOCITIES FOR STARTUP RUN");
1987 EETYPE("gen-vel", opts->bGenVel, yesno_names);
1988 RTYPE ("gen-temp", opts->tempi, 300.0);
1989 ITYPE ("gen-seed", opts->seed, -1);
1992 CCTYPE ("OPTIONS FOR BONDS");
1993 EETYPE("constraints", opts->nshake, constraints);
1994 CTYPE ("Type of constraint algorithm");
1995 EETYPE("constraint-algorithm", ir->eConstrAlg, econstr_names);
1996 CTYPE ("Do not constrain the start configuration");
1997 EETYPE("continuation", ir->bContinuation, yesno_names);
1998 CTYPE ("Use successive overrelaxation to reduce the number of shake iterations");
1999 EETYPE("Shake-SOR", ir->bShakeSOR, yesno_names);
2000 CTYPE ("Relative tolerance of shake");
2001 RTYPE ("shake-tol", ir->shake_tol, 0.0001);
2002 CTYPE ("Highest order in the expansion of the constraint coupling matrix");
2003 ITYPE ("lincs-order", ir->nProjOrder, 4);
2004 CTYPE ("Number of iterations in the final step of LINCS. 1 is fine for");
2005 CTYPE ("normal simulations, but use 2 to conserve energy in NVE runs.");
2006 CTYPE ("For energy minimization with constraints it should be 4 to 8.");
2007 ITYPE ("lincs-iter", ir->nLincsIter, 1);
2008 CTYPE ("Lincs will write a warning to the stderr if in one step a bond");
2009 CTYPE ("rotates over more degrees than");
2010 RTYPE ("lincs-warnangle", ir->LincsWarnAngle, 30.0);
2011 CTYPE ("Convert harmonic bonds to morse potentials");
2012 EETYPE("morse", opts->bMorse, yesno_names);
2014 /* Energy group exclusions */
2015 CCTYPE ("ENERGY GROUP EXCLUSIONS");
2016 CTYPE ("Pairs of energy groups for which all non-bonded interactions are excluded");
2017 STYPE ("energygrp-excl", is->egpexcl, NULL);
2021 CTYPE ("Number of walls, type, atom types, densities and box-z scale factor for Ewald");
2022 ITYPE ("nwall", ir->nwall, 0);
2023 EETYPE("wall-type", ir->wall_type, ewt_names);
2024 RTYPE ("wall-r-linpot", ir->wall_r_linpot, -1);
2025 STYPE ("wall-atomtype", is->wall_atomtype, NULL);
2026 STYPE ("wall-density", is->wall_density, NULL);
2027 RTYPE ("wall-ewald-zfac", ir->wall_ewald_zfac, 3);
2030 CCTYPE("COM PULLING");
2031 CTYPE("Pull type: no, umbrella, constraint or constant-force");
2032 EETYPE("pull", ir->ePull, epull_names);
2033 if (ir->ePull != epullNO)
2036 is->pull_grp = read_pullparams(&ninp, &inp, ir->pull, &opts->pull_start, wi);
2039 /* Enforced rotation */
2040 CCTYPE("ENFORCED ROTATION");
2041 CTYPE("Enforced rotation: No or Yes");
2042 EETYPE("rotation", ir->bRot, yesno_names);
2046 is->rot_grp = read_rotparams(&ninp, &inp, ir->rot, wi);
2050 CCTYPE("NMR refinement stuff");
2051 CTYPE ("Distance restraints type: No, Simple or Ensemble");
2052 EETYPE("disre", ir->eDisre, edisre_names);
2053 CTYPE ("Force weighting of pairs in one distance restraint: Conservative or Equal");
2054 EETYPE("disre-weighting", ir->eDisreWeighting, edisreweighting_names);
2055 CTYPE ("Use sqrt of the time averaged times the instantaneous violation");
2056 EETYPE("disre-mixed", ir->bDisreMixed, yesno_names);
2057 RTYPE ("disre-fc", ir->dr_fc, 1000.0);
2058 RTYPE ("disre-tau", ir->dr_tau, 0.0);
2059 CTYPE ("Output frequency for pair distances to energy file");
2060 ITYPE ("nstdisreout", ir->nstdisreout, 100);
2061 CTYPE ("Orientation restraints: No or Yes");
2062 EETYPE("orire", opts->bOrire, yesno_names);
2063 CTYPE ("Orientation restraints force constant and tau for time averaging");
2064 RTYPE ("orire-fc", ir->orires_fc, 0.0);
2065 RTYPE ("orire-tau", ir->orires_tau, 0.0);
2066 STYPE ("orire-fitgrp", is->orirefitgrp, NULL);
2067 CTYPE ("Output frequency for trace(SD) and S to energy file");
2068 ITYPE ("nstorireout", ir->nstorireout, 100);
2070 /* free energy variables */
2071 CCTYPE ("Free energy variables");
2072 EETYPE("free-energy", ir->efep, efep_names);
2073 STYPE ("couple-moltype", is->couple_moltype, NULL);
2074 EETYPE("couple-lambda0", opts->couple_lam0, couple_lam);
2075 EETYPE("couple-lambda1", opts->couple_lam1, couple_lam);
2076 EETYPE("couple-intramol", opts->bCoupleIntra, yesno_names);
2078 RTYPE ("init-lambda", fep->init_lambda, -1); /* start with -1 so
2080 it was not entered */
2081 ITYPE ("init-lambda-state", fep->init_fep_state, -1);
2082 RTYPE ("delta-lambda", fep->delta_lambda, 0.0);
2083 ITYPE ("nstdhdl", fep->nstdhdl, 50);
2084 STYPE ("fep-lambdas", is->fep_lambda[efptFEP], NULL);
2085 STYPE ("mass-lambdas", is->fep_lambda[efptMASS], NULL);
2086 STYPE ("coul-lambdas", is->fep_lambda[efptCOUL], NULL);
2087 STYPE ("vdw-lambdas", is->fep_lambda[efptVDW], NULL);
2088 STYPE ("bonded-lambdas", is->fep_lambda[efptBONDED], NULL);
2089 STYPE ("restraint-lambdas", is->fep_lambda[efptRESTRAINT], NULL);
2090 STYPE ("temperature-lambdas", is->fep_lambda[efptTEMPERATURE], NULL);
2091 ITYPE ("calc-lambda-neighbors", fep->lambda_neighbors, 1);
2092 STYPE ("init-lambda-weights", is->lambda_weights, NULL);
2093 EETYPE("dhdl-print-energy", fep->bPrintEnergy, yesno_names);
2094 RTYPE ("sc-alpha", fep->sc_alpha, 0.0);
2095 ITYPE ("sc-power", fep->sc_power, 1);
2096 RTYPE ("sc-r-power", fep->sc_r_power, 6.0);
2097 RTYPE ("sc-sigma", fep->sc_sigma, 0.3);
2098 EETYPE("sc-coul", fep->bScCoul, yesno_names);
2099 ITYPE ("dh_hist_size", fep->dh_hist_size, 0);
2100 RTYPE ("dh_hist_spacing", fep->dh_hist_spacing, 0.1);
2101 EETYPE("separate-dhdl-file", fep->separate_dhdl_file,
2102 separate_dhdl_file_names);
2103 EETYPE("dhdl-derivatives", fep->dhdl_derivatives, dhdl_derivatives_names);
2104 ITYPE ("dh_hist_size", fep->dh_hist_size, 0);
2105 RTYPE ("dh_hist_spacing", fep->dh_hist_spacing, 0.1);
2107 /* Non-equilibrium MD stuff */
2108 CCTYPE("Non-equilibrium MD stuff");
2109 STYPE ("acc-grps", is->accgrps, NULL);
2110 STYPE ("accelerate", is->acc, NULL);
2111 STYPE ("freezegrps", is->freeze, NULL);
2112 STYPE ("freezedim", is->frdim, NULL);
2113 RTYPE ("cos-acceleration", ir->cos_accel, 0);
2114 STYPE ("deform", is->deform, NULL);
2116 /* simulated tempering variables */
2117 CCTYPE("simulated tempering variables");
2118 EETYPE("simulated-tempering", ir->bSimTemp, yesno_names);
2119 EETYPE("simulated-tempering-scaling", ir->simtempvals->eSimTempScale, esimtemp_names);
2120 RTYPE("sim-temp-low", ir->simtempvals->simtemp_low, 300.0);
2121 RTYPE("sim-temp-high", ir->simtempvals->simtemp_high, 300.0);
2123 /* expanded ensemble variables */
2124 if (ir->efep == efepEXPANDED || ir->bSimTemp)
2126 read_expandedparams(&ninp, &inp, expand, wi);
2129 /* Electric fields */
2130 CCTYPE("Electric fields");
2131 CTYPE ("Format is number of terms (int) and for all terms an amplitude (real)");
2132 CTYPE ("and a phase angle (real)");
2133 STYPE ("E-x", is->efield_x, NULL);
2134 STYPE ("E-xt", is->efield_xt, NULL);
2135 STYPE ("E-y", is->efield_y, NULL);
2136 STYPE ("E-yt", is->efield_yt, NULL);
2137 STYPE ("E-z", is->efield_z, NULL);
2138 STYPE ("E-zt", is->efield_zt, NULL);
2140 CCTYPE("Ion/water position swapping for computational electrophysiology setups");
2141 CTYPE("Swap positions along direction: no, X, Y, Z");
2142 EETYPE("swapcoords", ir->eSwapCoords, eSwapTypes_names);
2143 if (ir->eSwapCoords != eswapNO)
2146 CTYPE("Swap attempt frequency");
2147 ITYPE("swap-frequency", ir->swap->nstswap, 1);
2148 CTYPE("Two index groups that contain the compartment-partitioning atoms");
2149 STYPE("split-group0", splitgrp0, NULL);
2150 STYPE("split-group1", splitgrp1, NULL);
2151 CTYPE("Use center of mass of split groups (yes/no), otherwise center of geometry is used");
2152 EETYPE("massw-split0", ir->swap->massw_split[0], yesno_names);
2153 EETYPE("massw-split1", ir->swap->massw_split[1], yesno_names);
2155 CTYPE("Group name of ions that can be exchanged with solvent molecules");
2156 STYPE("swap-group", swapgrp, NULL);
2157 CTYPE("Group name of solvent molecules");
2158 STYPE("solvent-group", solgrp, NULL);
2160 CTYPE("Split cylinder: radius, upper and lower extension (nm) (this will define the channels)");
2161 CTYPE("Note that the split cylinder settings do not have an influence on the swapping protocol,");
2162 CTYPE("however, if correctly defined, the ion permeation events are counted per channel");
2163 RTYPE("cyl0-r", ir->swap->cyl0r, 2.0);
2164 RTYPE("cyl0-up", ir->swap->cyl0u, 1.0);
2165 RTYPE("cyl0-down", ir->swap->cyl0l, 1.0);
2166 RTYPE("cyl1-r", ir->swap->cyl1r, 2.0);
2167 RTYPE("cyl1-up", ir->swap->cyl1u, 1.0);
2168 RTYPE("cyl1-down", ir->swap->cyl1l, 1.0);
2170 CTYPE("Average the number of ions per compartment over these many swap attempt steps");
2171 ITYPE("coupl-steps", ir->swap->nAverage, 10);
2172 CTYPE("Requested number of anions and cations for each of the two compartments");
2173 CTYPE("-1 means fix the numbers as found in time step 0");
2174 ITYPE("anionsA", ir->swap->nanions[0], -1);
2175 ITYPE("cationsA", ir->swap->ncations[0], -1);
2176 ITYPE("anionsB", ir->swap->nanions[1], -1);
2177 ITYPE("cationsB", ir->swap->ncations[1], -1);
2178 CTYPE("Start to swap ions if threshold difference to requested count is reached");
2179 RTYPE("threshold", ir->swap->threshold, 1.0);
2182 /* AdResS defined thingies */
2183 CCTYPE ("AdResS parameters");
2184 EETYPE("adress", ir->bAdress, yesno_names);
2187 snew(ir->adress, 1);
2188 read_adressparams(&ninp, &inp, ir->adress, wi);
2191 /* User defined thingies */
2192 CCTYPE ("User defined thingies");
2193 STYPE ("user1-grps", is->user1, NULL);
2194 STYPE ("user2-grps", is->user2, NULL);
2195 ITYPE ("userint1", ir->userint1, 0);
2196 ITYPE ("userint2", ir->userint2, 0);
2197 ITYPE ("userint3", ir->userint3, 0);
2198 ITYPE ("userint4", ir->userint4, 0);
2199 RTYPE ("userreal1", ir->userreal1, 0);
2200 RTYPE ("userreal2", ir->userreal2, 0);
2201 RTYPE ("userreal3", ir->userreal3, 0);
2202 RTYPE ("userreal4", ir->userreal4, 0);
2205 write_inpfile(mdparout, ninp, inp, FALSE, wi);
2206 for (i = 0; (i < ninp); i++)
2209 sfree(inp[i].value);
2213 /* Process options if necessary */
2214 for (m = 0; m < 2; m++)
2216 for (i = 0; i < 2*DIM; i++)
2225 if (sscanf(dumstr[m], "%lf", &(dumdub[m][XX])) != 1)
2227 warning_error(wi, "Pressure coupling not enough values (I need 1)");
2229 dumdub[m][YY] = dumdub[m][ZZ] = dumdub[m][XX];
2231 case epctSEMIISOTROPIC:
2232 case epctSURFACETENSION:
2233 if (sscanf(dumstr[m], "%lf%lf",
2234 &(dumdub[m][XX]), &(dumdub[m][ZZ])) != 2)
2236 warning_error(wi, "Pressure coupling not enough values (I need 2)");
2238 dumdub[m][YY] = dumdub[m][XX];
2240 case epctANISOTROPIC:
2241 if (sscanf(dumstr[m], "%lf%lf%lf%lf%lf%lf",
2242 &(dumdub[m][XX]), &(dumdub[m][YY]), &(dumdub[m][ZZ]),
2243 &(dumdub[m][3]), &(dumdub[m][4]), &(dumdub[m][5])) != 6)
2245 warning_error(wi, "Pressure coupling not enough values (I need 6)");
2249 gmx_fatal(FARGS, "Pressure coupling type %s not implemented yet",
2250 epcoupltype_names[ir->epct]);
2254 clear_mat(ir->ref_p);
2255 clear_mat(ir->compress);
2256 for (i = 0; i < DIM; i++)
2258 ir->ref_p[i][i] = dumdub[1][i];
2259 ir->compress[i][i] = dumdub[0][i];
2261 if (ir->epct == epctANISOTROPIC)
2263 ir->ref_p[XX][YY] = dumdub[1][3];
2264 ir->ref_p[XX][ZZ] = dumdub[1][4];
2265 ir->ref_p[YY][ZZ] = dumdub[1][5];
2266 if (ir->ref_p[XX][YY] != 0 && ir->ref_p[XX][ZZ] != 0 && ir->ref_p[YY][ZZ] != 0)
2268 warning(wi, "All off-diagonal reference pressures are non-zero. Are you sure you want to apply a threefold shear stress?\n");
2270 ir->compress[XX][YY] = dumdub[0][3];
2271 ir->compress[XX][ZZ] = dumdub[0][4];
2272 ir->compress[YY][ZZ] = dumdub[0][5];
2273 for (i = 0; i < DIM; i++)
2275 for (m = 0; m < i; m++)
2277 ir->ref_p[i][m] = ir->ref_p[m][i];
2278 ir->compress[i][m] = ir->compress[m][i];
2283 if (ir->comm_mode == ecmNO)
2288 opts->couple_moltype = NULL;
2289 if (strlen(is->couple_moltype) > 0)
2291 if (ir->efep != efepNO)
2293 opts->couple_moltype = strdup(is->couple_moltype);
2294 if (opts->couple_lam0 == opts->couple_lam1)
2296 warning(wi, "The lambda=0 and lambda=1 states for coupling are identical");
2298 if (ir->eI == eiMD && (opts->couple_lam0 == ecouplamNONE ||
2299 opts->couple_lam1 == ecouplamNONE))
2301 warning(wi, "For proper sampling of the (nearly) decoupled state, stochastic dynamics should be used");
2306 warning(wi, "Can not couple a molecule with free_energy = no");
2309 /* FREE ENERGY AND EXPANDED ENSEMBLE OPTIONS */
2310 if (ir->efep != efepNO)
2312 if (fep->delta_lambda > 0)
2314 ir->efep = efepSLOWGROWTH;
2320 fep->bPrintEnergy = TRUE;
2321 /* always print out the energy to dhdl if we are doing expanded ensemble, since we need the total energy
2322 if the temperature is changing. */
2325 if ((ir->efep != efepNO) || ir->bSimTemp)
2327 ir->bExpanded = FALSE;
2328 if ((ir->efep == efepEXPANDED) || ir->bSimTemp)
2330 ir->bExpanded = TRUE;
2332 do_fep_params(ir, is->fep_lambda, is->lambda_weights);
2333 if (ir->bSimTemp) /* done after fep params */
2335 do_simtemp_params(ir);
2340 ir->fepvals->n_lambda = 0;
2343 /* WALL PARAMETERS */
2345 do_wall_params(ir, is->wall_atomtype, is->wall_density, opts);
2347 /* ORIENTATION RESTRAINT PARAMETERS */
2349 if (opts->bOrire && str_nelem(is->orirefitgrp, MAXPTR, NULL) != 1)
2351 warning_error(wi, "ERROR: Need one orientation restraint fit group\n");
2354 /* DEFORMATION PARAMETERS */
2356 clear_mat(ir->deform);
2357 for (i = 0; i < 6; i++)
2361 m = sscanf(is->deform, "%lf %lf %lf %lf %lf %lf",
2362 &(dumdub[0][0]), &(dumdub[0][1]), &(dumdub[0][2]),
2363 &(dumdub[0][3]), &(dumdub[0][4]), &(dumdub[0][5]));
2364 for (i = 0; i < 3; i++)
2366 ir->deform[i][i] = dumdub[0][i];
2368 ir->deform[YY][XX] = dumdub[0][3];
2369 ir->deform[ZZ][XX] = dumdub[0][4];
2370 ir->deform[ZZ][YY] = dumdub[0][5];
2371 if (ir->epc != epcNO)
2373 for (i = 0; i < 3; i++)
2375 for (j = 0; j <= i; j++)
2377 if (ir->deform[i][j] != 0 && ir->compress[i][j] != 0)
2379 warning_error(wi, "A box element has deform set and compressibility > 0");
2383 for (i = 0; i < 3; i++)
2385 for (j = 0; j < i; j++)
2387 if (ir->deform[i][j] != 0)
2389 for (m = j; m < DIM; m++)
2391 if (ir->compress[m][j] != 0)
2393 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.");
2394 warning(wi, warn_buf);
2402 /* Ion/water position swapping checks */
2403 if (ir->eSwapCoords != eswapNO)
2405 if (ir->swap->nstswap < 1)
2407 warning_error(wi, "swap_frequency must be 1 or larger when ion swapping is requested");
2409 if (ir->swap->nAverage < 1)
2411 warning_error(wi, "coupl_steps must be 1 or larger.\n");
2413 if (ir->swap->threshold < 1.0)
2415 warning_error(wi, "Ion count threshold must be at least 1.\n");
2423 static int search_QMstring(const char *s, int ng, const char *gn[])
2425 /* same as normal search_string, but this one searches QM strings */
2428 for (i = 0; (i < ng); i++)
2430 if (gmx_strcasecmp(s, gn[i]) == 0)
2436 gmx_fatal(FARGS, "this QM method or basisset (%s) is not implemented\n!", s);
2440 } /* search_QMstring */
2442 /* We would like gn to be const as well, but C doesn't allow this */
2443 int search_string(const char *s, int ng, char *gn[])
2447 for (i = 0; (i < ng); i++)
2449 if (gmx_strcasecmp(s, gn[i]) == 0)
2456 "Group %s referenced in the .mdp file was not found in the index file.\n"
2457 "Group names must match either [moleculetype] names or custom index group\n"
2458 "names, in which case you must supply an index file to the '-n' option\n"
2465 static gmx_bool do_numbering(int natoms, gmx_groups_t *groups, int ng, char *ptrs[],
2466 t_blocka *block, char *gnames[],
2467 int gtype, int restnm,
2468 int grptp, gmx_bool bVerbose,
2471 unsigned short *cbuf;
2472 t_grps *grps = &(groups->grps[gtype]);
2473 int i, j, gid, aj, ognr, ntot = 0;
2476 char warn_buf[STRLEN];
2480 fprintf(debug, "Starting numbering %d groups of type %d\n", ng, gtype);
2483 title = gtypes[gtype];
2486 /* Mark all id's as not set */
2487 for (i = 0; (i < natoms); i++)
2492 snew(grps->nm_ind, ng+1); /* +1 for possible rest group */
2493 for (i = 0; (i < ng); i++)
2495 /* Lookup the group name in the block structure */
2496 gid = search_string(ptrs[i], block->nr, gnames);
2497 if ((grptp != egrptpONE) || (i == 0))
2499 grps->nm_ind[grps->nr++] = gid;
2503 fprintf(debug, "Found gid %d for group %s\n", gid, ptrs[i]);
2506 /* Now go over the atoms in the group */
2507 for (j = block->index[gid]; (j < block->index[gid+1]); j++)
2512 /* Range checking */
2513 if ((aj < 0) || (aj >= natoms))
2515 gmx_fatal(FARGS, "Invalid atom number %d in indexfile", aj);
2517 /* Lookup up the old group number */
2521 gmx_fatal(FARGS, "Atom %d in multiple %s groups (%d and %d)",
2522 aj+1, title, ognr+1, i+1);
2526 /* Store the group number in buffer */
2527 if (grptp == egrptpONE)
2540 /* Now check whether we have done all atoms */
2544 if (grptp == egrptpALL)
2546 gmx_fatal(FARGS, "%d atoms are not part of any of the %s groups",
2547 natoms-ntot, title);
2549 else if (grptp == egrptpPART)
2551 sprintf(warn_buf, "%d atoms are not part of any of the %s groups",
2552 natoms-ntot, title);
2553 warning_note(wi, warn_buf);
2555 /* Assign all atoms currently unassigned to a rest group */
2556 for (j = 0; (j < natoms); j++)
2558 if (cbuf[j] == NOGID)
2564 if (grptp != egrptpPART)
2569 "Making dummy/rest group for %s containing %d elements\n",
2570 title, natoms-ntot);
2572 /* Add group name "rest" */
2573 grps->nm_ind[grps->nr] = restnm;
2575 /* Assign the rest name to all atoms not currently assigned to a group */
2576 for (j = 0; (j < natoms); j++)
2578 if (cbuf[j] == NOGID)
2587 if (grps->nr == 1 && (ntot == 0 || ntot == natoms))
2589 /* All atoms are part of one (or no) group, no index required */
2590 groups->ngrpnr[gtype] = 0;
2591 groups->grpnr[gtype] = NULL;
2595 groups->ngrpnr[gtype] = natoms;
2596 snew(groups->grpnr[gtype], natoms);
2597 for (j = 0; (j < natoms); j++)
2599 groups->grpnr[gtype][j] = cbuf[j];
2605 return (bRest && grptp == egrptpPART);
2608 static void calc_nrdf(gmx_mtop_t *mtop, t_inputrec *ir, char **gnames)
2611 gmx_groups_t *groups;
2613 int natoms, ai, aj, i, j, d, g, imin, jmin;
2615 int *nrdf2, *na_vcm, na_tot;
2616 double *nrdf_tc, *nrdf_vcm, nrdf_uc, n_sub = 0;
2617 gmx_mtop_atomloop_all_t aloop;
2619 int mb, mol, ftype, as;
2620 gmx_molblock_t *molb;
2621 gmx_moltype_t *molt;
2624 * First calc 3xnr-atoms for each group
2625 * then subtract half a degree of freedom for each constraint
2627 * Only atoms and nuclei contribute to the degrees of freedom...
2632 groups = &mtop->groups;
2633 natoms = mtop->natoms;
2635 /* Allocate one more for a possible rest group */
2636 /* We need to sum degrees of freedom into doubles,
2637 * since floats give too low nrdf's above 3 million atoms.
2639 snew(nrdf_tc, groups->grps[egcTC].nr+1);
2640 snew(nrdf_vcm, groups->grps[egcVCM].nr+1);
2641 snew(na_vcm, groups->grps[egcVCM].nr+1);
2643 for (i = 0; i < groups->grps[egcTC].nr; i++)
2647 for (i = 0; i < groups->grps[egcVCM].nr+1; i++)
2652 snew(nrdf2, natoms);
2653 aloop = gmx_mtop_atomloop_all_init(mtop);
2654 while (gmx_mtop_atomloop_all_next(aloop, &i, &atom))
2657 if (atom->ptype == eptAtom || atom->ptype == eptNucleus)
2659 g = ggrpnr(groups, egcFREEZE, i);
2660 /* Double count nrdf for particle i */
2661 for (d = 0; d < DIM; d++)
2663 if (opts->nFreeze[g][d] == 0)
2668 nrdf_tc [ggrpnr(groups, egcTC, i)] += 0.5*nrdf2[i];
2669 nrdf_vcm[ggrpnr(groups, egcVCM, i)] += 0.5*nrdf2[i];
2674 for (mb = 0; mb < mtop->nmolblock; mb++)
2676 molb = &mtop->molblock[mb];
2677 molt = &mtop->moltype[molb->type];
2678 atom = molt->atoms.atom;
2679 for (mol = 0; mol < molb->nmol; mol++)
2681 for (ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++)
2683 ia = molt->ilist[ftype].iatoms;
2684 for (i = 0; i < molt->ilist[ftype].nr; )
2686 /* Subtract degrees of freedom for the constraints,
2687 * if the particles still have degrees of freedom left.
2688 * If one of the particles is a vsite or a shell, then all
2689 * constraint motion will go there, but since they do not
2690 * contribute to the constraints the degrees of freedom do not
2695 if (((atom[ia[1]].ptype == eptNucleus) ||
2696 (atom[ia[1]].ptype == eptAtom)) &&
2697 ((atom[ia[2]].ptype == eptNucleus) ||
2698 (atom[ia[2]].ptype == eptAtom)))
2716 imin = min(imin, nrdf2[ai]);
2717 jmin = min(jmin, nrdf2[aj]);
2720 nrdf_tc [ggrpnr(groups, egcTC, ai)] -= 0.5*imin;
2721 nrdf_tc [ggrpnr(groups, egcTC, aj)] -= 0.5*jmin;
2722 nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5*imin;
2723 nrdf_vcm[ggrpnr(groups, egcVCM, aj)] -= 0.5*jmin;
2725 ia += interaction_function[ftype].nratoms+1;
2726 i += interaction_function[ftype].nratoms+1;
2729 ia = molt->ilist[F_SETTLE].iatoms;
2730 for (i = 0; i < molt->ilist[F_SETTLE].nr; )
2732 /* Subtract 1 dof from every atom in the SETTLE */
2733 for (j = 0; j < 3; j++)
2736 imin = min(2, nrdf2[ai]);
2738 nrdf_tc [ggrpnr(groups, egcTC, ai)] -= 0.5*imin;
2739 nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5*imin;
2744 as += molt->atoms.nr;
2748 if (ir->ePull == epullCONSTRAINT)
2750 /* Correct nrdf for the COM constraints.
2751 * We correct using the TC and VCM group of the first atom
2752 * in the reference and pull group. If atoms in one pull group
2753 * belong to different TC or VCM groups it is anyhow difficult
2754 * to determine the optimal nrdf assignment.
2758 for (i = 0; i < pull->ncoord; i++)
2762 for (j = 0; j < 2; j++)
2764 const t_pull_group *pgrp;
2766 pgrp = &pull->group[pull->coord[i].group[j]];
2770 /* Subtract 1/2 dof from each group */
2772 nrdf_tc [ggrpnr(groups, egcTC, ai)] -= 0.5*imin;
2773 nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5*imin;
2774 if (nrdf_tc[ggrpnr(groups, egcTC, ai)] < 0)
2776 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)]]);
2781 /* We need to subtract the whole DOF from group j=1 */
2788 if (ir->nstcomm != 0)
2790 /* Subtract 3 from the number of degrees of freedom in each vcm group
2791 * when com translation is removed and 6 when rotation is removed
2794 switch (ir->comm_mode)
2797 n_sub = ndof_com(ir);
2804 gmx_incons("Checking comm_mode");
2807 for (i = 0; i < groups->grps[egcTC].nr; i++)
2809 /* Count the number of atoms of TC group i for every VCM group */
2810 for (j = 0; j < groups->grps[egcVCM].nr+1; j++)
2815 for (ai = 0; ai < natoms; ai++)
2817 if (ggrpnr(groups, egcTC, ai) == i)
2819 na_vcm[ggrpnr(groups, egcVCM, ai)]++;
2823 /* Correct for VCM removal according to the fraction of each VCM
2824 * group present in this TC group.
2826 nrdf_uc = nrdf_tc[i];
2829 fprintf(debug, "T-group[%d] nrdf_uc = %g, n_sub = %g\n",
2833 for (j = 0; j < groups->grps[egcVCM].nr+1; j++)
2835 if (nrdf_vcm[j] > n_sub)
2837 nrdf_tc[i] += nrdf_uc*((double)na_vcm[j]/(double)na_tot)*
2838 (nrdf_vcm[j] - n_sub)/nrdf_vcm[j];
2842 fprintf(debug, " nrdf_vcm[%d] = %g, nrdf = %g\n",
2843 j, nrdf_vcm[j], nrdf_tc[i]);
2848 for (i = 0; (i < groups->grps[egcTC].nr); i++)
2850 opts->nrdf[i] = nrdf_tc[i];
2851 if (opts->nrdf[i] < 0)
2856 "Number of degrees of freedom in T-Coupling group %s is %.2f\n",
2857 gnames[groups->grps[egcTC].nm_ind[i]], opts->nrdf[i]);
2866 static void decode_cos(char *s, t_cosines *cosine)
2869 char format[STRLEN], f1[STRLEN];
2881 sscanf(t, "%d", &(cosine->n));
2888 snew(cosine->a, cosine->n);
2889 snew(cosine->phi, cosine->n);
2891 sprintf(format, "%%*d");
2892 for (i = 0; (i < cosine->n); i++)
2895 strcat(f1, "%lf%lf");
2896 if (sscanf(t, f1, &a, &phi) < 2)
2898 gmx_fatal(FARGS, "Invalid input for electric field shift: '%s'", t);
2901 cosine->phi[i] = phi;
2902 strcat(format, "%*lf%*lf");
2909 static gmx_bool do_egp_flag(t_inputrec *ir, gmx_groups_t *groups,
2910 const char *option, const char *val, int flag)
2912 /* The maximum number of energy group pairs would be MAXPTR*(MAXPTR+1)/2.
2913 * But since this is much larger than STRLEN, such a line can not be parsed.
2914 * The real maximum is the number of names that fit in a string: STRLEN/2.
2916 #define EGP_MAX (STRLEN/2)
2917 int nelem, i, j, k, nr;
2918 char *names[EGP_MAX];
2922 gnames = groups->grpname;
2924 nelem = str_nelem(val, EGP_MAX, names);
2927 gmx_fatal(FARGS, "The number of groups for %s is odd", option);
2929 nr = groups->grps[egcENER].nr;
2931 for (i = 0; i < nelem/2; i++)
2935 gmx_strcasecmp(names[2*i], *(gnames[groups->grps[egcENER].nm_ind[j]])))
2941 gmx_fatal(FARGS, "%s in %s is not an energy group\n",
2942 names[2*i], option);
2946 gmx_strcasecmp(names[2*i+1], *(gnames[groups->grps[egcENER].nm_ind[k]])))
2952 gmx_fatal(FARGS, "%s in %s is not an energy group\n",
2953 names[2*i+1], option);
2955 if ((j < nr) && (k < nr))
2957 ir->opts.egp_flags[nr*j+k] |= flag;
2958 ir->opts.egp_flags[nr*k+j] |= flag;
2967 static void make_swap_groups(
2976 int ig = -1, i = 0, j;
2980 /* Just a quick check here, more thorough checks are in mdrun */
2981 if (strcmp(splitg0name, splitg1name) == 0)
2983 gmx_fatal(FARGS, "The split groups can not both be '%s'.", splitg0name);
2986 /* First get the swap group index atoms */
2987 ig = search_string(swapgname, grps->nr, gnames);
2988 swap->nat = grps->index[ig+1] - grps->index[ig];
2991 fprintf(stderr, "Swap group '%s' contains %d atoms.\n", swapgname, swap->nat);
2992 snew(swap->ind, swap->nat);
2993 for (i = 0; i < swap->nat; i++)
2995 swap->ind[i] = grps->a[grps->index[ig]+i];
3000 gmx_fatal(FARGS, "You defined an empty group of atoms for swapping.");
3003 /* Now do so for the split groups */
3004 for (j = 0; j < 2; j++)
3008 splitg = splitg0name;
3012 splitg = splitg1name;
3015 ig = search_string(splitg, grps->nr, gnames);
3016 swap->nat_split[j] = grps->index[ig+1] - grps->index[ig];
3017 if (swap->nat_split[j] > 0)
3019 fprintf(stderr, "Split group %d '%s' contains %d atom%s.\n",
3020 j, splitg, swap->nat_split[j], (swap->nat_split[j] > 1) ? "s" : "");
3021 snew(swap->ind_split[j], swap->nat_split[j]);
3022 for (i = 0; i < swap->nat_split[j]; i++)
3024 swap->ind_split[j][i] = grps->a[grps->index[ig]+i];
3029 gmx_fatal(FARGS, "Split group %d has to contain at least 1 atom!", j);
3033 /* Now get the solvent group index atoms */
3034 ig = search_string(solgname, grps->nr, gnames);
3035 swap->nat_sol = grps->index[ig+1] - grps->index[ig];
3036 if (swap->nat_sol > 0)
3038 fprintf(stderr, "Solvent group '%s' contains %d atoms.\n", solgname, swap->nat_sol);
3039 snew(swap->ind_sol, swap->nat_sol);
3040 for (i = 0; i < swap->nat_sol; i++)
3042 swap->ind_sol[i] = grps->a[grps->index[ig]+i];
3047 gmx_fatal(FARGS, "You defined an empty group of solvent. Cannot exchange ions.");
3052 void do_index(const char* mdparin, const char *ndx,
3055 t_inputrec *ir, rvec *v,
3059 gmx_groups_t *groups;
3063 char warnbuf[STRLEN], **gnames;
3064 int nr, ntcg, ntau_t, nref_t, nacc, nofg, nSA, nSA_points, nSA_time, nSA_temp;
3067 int nacg, nfreeze, nfrdim, nenergy, nvcm, nuser;
3068 char *ptr1[MAXPTR], *ptr2[MAXPTR], *ptr3[MAXPTR];
3069 int i, j, k, restnm;
3071 gmx_bool bExcl, bTable, bSetTCpar, bAnneal, bRest;
3072 int nQMmethod, nQMbasis, nQMcharge, nQMmult, nbSH, nCASorb, nCASelec,
3073 nSAon, nSAoff, nSAsteps, nQMg, nbOPT, nbTS;
3074 char warn_buf[STRLEN];
3078 fprintf(stderr, "processing index file...\n");
3084 snew(grps->index, 1);
3086 atoms_all = gmx_mtop_global_atoms(mtop);
3087 analyse(&atoms_all, grps, &gnames, FALSE, TRUE);
3088 free_t_atoms(&atoms_all, FALSE);
3092 grps = init_index(ndx, &gnames);
3095 groups = &mtop->groups;
3096 natoms = mtop->natoms;
3097 symtab = &mtop->symtab;
3099 snew(groups->grpname, grps->nr+1);
3101 for (i = 0; (i < grps->nr); i++)
3103 groups->grpname[i] = put_symtab(symtab, gnames[i]);
3105 groups->grpname[i] = put_symtab(symtab, "rest");
3107 srenew(gnames, grps->nr+1);
3108 gnames[restnm] = *(groups->grpname[i]);
3109 groups->ngrpname = grps->nr+1;
3111 set_warning_line(wi, mdparin, -1);
3113 ntau_t = str_nelem(is->tau_t, MAXPTR, ptr1);
3114 nref_t = str_nelem(is->ref_t, MAXPTR, ptr2);
3115 ntcg = str_nelem(is->tcgrps, MAXPTR, ptr3);
3116 if ((ntau_t != ntcg) || (nref_t != ntcg))
3118 gmx_fatal(FARGS, "Invalid T coupling input: %d groups, %d ref-t values and "
3119 "%d tau-t values", ntcg, nref_t, ntau_t);
3122 bSetTCpar = (ir->etc || EI_SD(ir->eI) || ir->eI == eiBD || EI_TPI(ir->eI));
3123 do_numbering(natoms, groups, ntcg, ptr3, grps, gnames, egcTC,
3124 restnm, bSetTCpar ? egrptpALL : egrptpALL_GENREST, bVerbose, wi);
3125 nr = groups->grps[egcTC].nr;
3127 snew(ir->opts.nrdf, nr);
3128 snew(ir->opts.tau_t, nr);
3129 snew(ir->opts.ref_t, nr);
3130 if (ir->eI == eiBD && ir->bd_fric == 0)
3132 fprintf(stderr, "bd-fric=0, so tau-t will be used as the inverse friction constant(s)\n");
3139 gmx_fatal(FARGS, "Not enough ref-t and tau-t values!");
3143 for (i = 0; (i < nr); i++)
3145 ir->opts.tau_t[i] = strtod(ptr1[i], NULL);
3146 if ((ir->eI == eiBD || ir->eI == eiSD2) && ir->opts.tau_t[i] <= 0)
3148 sprintf(warn_buf, "With integrator %s tau-t should be larger than 0", ei_names[ir->eI]);
3149 warning_error(wi, warn_buf);
3152 if (ir->etc != etcVRESCALE && ir->opts.tau_t[i] == 0)
3154 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.");
3157 if (ir->opts.tau_t[i] >= 0)
3159 tau_min = min(tau_min, ir->opts.tau_t[i]);
3162 if (ir->etc != etcNO && ir->nsttcouple == -1)
3164 ir->nsttcouple = ir_optimal_nsttcouple(ir);
3169 if ((ir->etc == etcNOSEHOOVER) && (ir->epc == epcBERENDSEN))
3171 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");
3173 if ((ir->epc == epcMTTK) && (ir->etc > etcNO))
3175 if (ir->nstpcouple != ir->nsttcouple)
3177 int mincouple = min(ir->nstpcouple, ir->nsttcouple);
3178 ir->nstpcouple = ir->nsttcouple = mincouple;
3179 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);
3180 warning_note(wi, warn_buf);
3184 /* velocity verlet with averaged kinetic energy KE = 0.5*(v(t+1/2) - v(t-1/2)) is implemented
3185 primarily for testing purposes, and does not work with temperature coupling other than 1 */
3187 if (ETC_ANDERSEN(ir->etc))
3189 if (ir->nsttcouple != 1)
3192 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");
3193 warning_note(wi, warn_buf);
3196 nstcmin = tcouple_min_integration_steps(ir->etc);
3199 if (tau_min/(ir->delta_t*ir->nsttcouple) < nstcmin)
3201 sprintf(warn_buf, "For proper integration of the %s thermostat, tau-t (%g) should be at least %d times larger than nsttcouple*dt (%g)",
3202 ETCOUPLTYPE(ir->etc),
3204 ir->nsttcouple*ir->delta_t);
3205 warning(wi, warn_buf);
3208 for (i = 0; (i < nr); i++)
3210 ir->opts.ref_t[i] = strtod(ptr2[i], NULL);
3211 if (ir->opts.ref_t[i] < 0)
3213 gmx_fatal(FARGS, "ref-t for group %d negative", i);
3216 /* set the lambda mc temperature to the md integrator temperature (which should be defined
3217 if we are in this conditional) if mc_temp is negative */
3218 if (ir->expandedvals->mc_temp < 0)
3220 ir->expandedvals->mc_temp = ir->opts.ref_t[0]; /*for now, set to the first reft */
3224 /* Simulated annealing for each group. There are nr groups */
3225 nSA = str_nelem(is->anneal, MAXPTR, ptr1);
3226 if (nSA == 1 && (ptr1[0][0] == 'n' || ptr1[0][0] == 'N'))
3230 if (nSA > 0 && nSA != nr)
3232 gmx_fatal(FARGS, "Not enough annealing values: %d (for %d groups)\n", nSA, nr);
3236 snew(ir->opts.annealing, nr);
3237 snew(ir->opts.anneal_npoints, nr);
3238 snew(ir->opts.anneal_time, nr);
3239 snew(ir->opts.anneal_temp, nr);
3240 for (i = 0; i < nr; i++)
3242 ir->opts.annealing[i] = eannNO;
3243 ir->opts.anneal_npoints[i] = 0;
3244 ir->opts.anneal_time[i] = NULL;
3245 ir->opts.anneal_temp[i] = NULL;
3250 for (i = 0; i < nr; i++)
3252 if (ptr1[i][0] == 'n' || ptr1[i][0] == 'N')
3254 ir->opts.annealing[i] = eannNO;
3256 else if (ptr1[i][0] == 's' || ptr1[i][0] == 'S')
3258 ir->opts.annealing[i] = eannSINGLE;
3261 else if (ptr1[i][0] == 'p' || ptr1[i][0] == 'P')
3263 ir->opts.annealing[i] = eannPERIODIC;
3269 /* Read the other fields too */
3270 nSA_points = str_nelem(is->anneal_npoints, MAXPTR, ptr1);
3271 if (nSA_points != nSA)
3273 gmx_fatal(FARGS, "Found %d annealing-npoints values for %d groups\n", nSA_points, nSA);
3275 for (k = 0, i = 0; i < nr; i++)
3277 ir->opts.anneal_npoints[i] = strtol(ptr1[i], NULL, 10);
3278 if (ir->opts.anneal_npoints[i] == 1)
3280 gmx_fatal(FARGS, "Please specify at least a start and an end point for annealing\n");
3282 snew(ir->opts.anneal_time[i], ir->opts.anneal_npoints[i]);
3283 snew(ir->opts.anneal_temp[i], ir->opts.anneal_npoints[i]);
3284 k += ir->opts.anneal_npoints[i];
3287 nSA_time = str_nelem(is->anneal_time, MAXPTR, ptr1);
3290 gmx_fatal(FARGS, "Found %d annealing-time values, wanter %d\n", nSA_time, k);
3292 nSA_temp = str_nelem(is->anneal_temp, MAXPTR, ptr2);
3295 gmx_fatal(FARGS, "Found %d annealing-temp values, wanted %d\n", nSA_temp, k);
3298 for (i = 0, k = 0; i < nr; i++)
3301 for (j = 0; j < ir->opts.anneal_npoints[i]; j++)
3303 ir->opts.anneal_time[i][j] = strtod(ptr1[k], NULL);
3304 ir->opts.anneal_temp[i][j] = strtod(ptr2[k], NULL);
3307 if (ir->opts.anneal_time[i][0] > (ir->init_t+GMX_REAL_EPS))
3309 gmx_fatal(FARGS, "First time point for annealing > init_t.\n");
3315 if (ir->opts.anneal_time[i][j] < ir->opts.anneal_time[i][j-1])
3317 gmx_fatal(FARGS, "Annealing timepoints out of order: t=%f comes after t=%f\n",
3318 ir->opts.anneal_time[i][j], ir->opts.anneal_time[i][j-1]);
3321 if (ir->opts.anneal_temp[i][j] < 0)
3323 gmx_fatal(FARGS, "Found negative temperature in annealing: %f\n", ir->opts.anneal_temp[i][j]);
3328 /* Print out some summary information, to make sure we got it right */
3329 for (i = 0, k = 0; i < nr; i++)
3331 if (ir->opts.annealing[i] != eannNO)
3333 j = groups->grps[egcTC].nm_ind[i];
3334 fprintf(stderr, "Simulated annealing for group %s: %s, %d timepoints\n",
3335 *(groups->grpname[j]), eann_names[ir->opts.annealing[i]],
3336 ir->opts.anneal_npoints[i]);
3337 fprintf(stderr, "Time (ps) Temperature (K)\n");
3338 /* All terms except the last one */
3339 for (j = 0; j < (ir->opts.anneal_npoints[i]-1); j++)
3341 fprintf(stderr, "%9.1f %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3344 /* Finally the last one */
3345 j = ir->opts.anneal_npoints[i]-1;
3346 if (ir->opts.annealing[i] == eannSINGLE)
3348 fprintf(stderr, "%9.1f- %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3352 fprintf(stderr, "%9.1f %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3353 if (fabs(ir->opts.anneal_temp[i][j]-ir->opts.anneal_temp[i][0]) > GMX_REAL_EPS)
3355 warning_note(wi, "There is a temperature jump when your annealing loops back.\n");
3364 if (ir->ePull != epullNO)
3366 make_pull_groups(ir->pull, is->pull_grp, grps, gnames);
3368 make_pull_coords(ir->pull);
3373 make_rotation_groups(ir->rot, is->rot_grp, grps, gnames);
3376 if (ir->eSwapCoords != eswapNO)
3378 make_swap_groups(ir->swap, swapgrp, splitgrp0, splitgrp1, solgrp, grps, gnames);
3381 nacc = str_nelem(is->acc, MAXPTR, ptr1);
3382 nacg = str_nelem(is->accgrps, MAXPTR, ptr2);
3383 if (nacg*DIM != nacc)
3385 gmx_fatal(FARGS, "Invalid Acceleration input: %d groups and %d acc. values",
3388 do_numbering(natoms, groups, nacg, ptr2, grps, gnames, egcACC,
3389 restnm, egrptpALL_GENREST, bVerbose, wi);
3390 nr = groups->grps[egcACC].nr;
3391 snew(ir->opts.acc, nr);
3392 ir->opts.ngacc = nr;
3394 for (i = k = 0; (i < nacg); i++)
3396 for (j = 0; (j < DIM); j++, k++)
3398 ir->opts.acc[i][j] = strtod(ptr1[k], NULL);
3401 for (; (i < nr); i++)
3403 for (j = 0; (j < DIM); j++)
3405 ir->opts.acc[i][j] = 0;
3409 nfrdim = str_nelem(is->frdim, MAXPTR, ptr1);
3410 nfreeze = str_nelem(is->freeze, MAXPTR, ptr2);
3411 if (nfrdim != DIM*nfreeze)
3413 gmx_fatal(FARGS, "Invalid Freezing input: %d groups and %d freeze values",
3416 do_numbering(natoms, groups, nfreeze, ptr2, grps, gnames, egcFREEZE,
3417 restnm, egrptpALL_GENREST, bVerbose, wi);
3418 nr = groups->grps[egcFREEZE].nr;
3419 ir->opts.ngfrz = nr;
3420 snew(ir->opts.nFreeze, nr);
3421 for (i = k = 0; (i < nfreeze); i++)
3423 for (j = 0; (j < DIM); j++, k++)
3425 ir->opts.nFreeze[i][j] = (gmx_strncasecmp(ptr1[k], "Y", 1) == 0);
3426 if (!ir->opts.nFreeze[i][j])
3428 if (gmx_strncasecmp(ptr1[k], "N", 1) != 0)
3430 sprintf(warnbuf, "Please use Y(ES) or N(O) for freezedim only "
3431 "(not %s)", ptr1[k]);
3432 warning(wi, warn_buf);
3437 for (; (i < nr); i++)
3439 for (j = 0; (j < DIM); j++)
3441 ir->opts.nFreeze[i][j] = 0;
3445 nenergy = str_nelem(is->energy, MAXPTR, ptr1);
3446 do_numbering(natoms, groups, nenergy, ptr1, grps, gnames, egcENER,
3447 restnm, egrptpALL_GENREST, bVerbose, wi);
3448 add_wall_energrps(groups, ir->nwall, symtab);
3449 ir->opts.ngener = groups->grps[egcENER].nr;
3450 nvcm = str_nelem(is->vcm, MAXPTR, ptr1);
3452 do_numbering(natoms, groups, nvcm, ptr1, grps, gnames, egcVCM,
3453 restnm, nvcm == 0 ? egrptpALL_GENREST : egrptpPART, bVerbose, wi);
3456 warning(wi, "Some atoms are not part of any center of mass motion removal group.\n"
3457 "This may lead to artifacts.\n"
3458 "In most cases one should use one group for the whole system.");
3461 /* Now we have filled the freeze struct, so we can calculate NRDF */
3462 calc_nrdf(mtop, ir, gnames);
3468 /* Must check per group! */
3469 for (i = 0; (i < ir->opts.ngtc); i++)
3471 ntot += ir->opts.nrdf[i];
3473 if (ntot != (DIM*natoms))
3475 fac = sqrt(ntot/(DIM*natoms));
3478 fprintf(stderr, "Scaling velocities by a factor of %.3f to account for constraints\n"
3479 "and removal of center of mass motion\n", fac);
3481 for (i = 0; (i < natoms); i++)
3483 svmul(fac, v[i], v[i]);
3488 nuser = str_nelem(is->user1, MAXPTR, ptr1);
3489 do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcUser1,
3490 restnm, egrptpALL_GENREST, bVerbose, wi);
3491 nuser = str_nelem(is->user2, MAXPTR, ptr1);
3492 do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcUser2,
3493 restnm, egrptpALL_GENREST, bVerbose, wi);
3494 nuser = str_nelem(is->x_compressed_groups, MAXPTR, ptr1);
3495 do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcCompressedX,
3496 restnm, egrptpONE, bVerbose, wi);
3497 nofg = str_nelem(is->orirefitgrp, MAXPTR, ptr1);
3498 do_numbering(natoms, groups, nofg, ptr1, grps, gnames, egcORFIT,
3499 restnm, egrptpALL_GENREST, bVerbose, wi);
3501 /* QMMM input processing */
3502 nQMg = str_nelem(is->QMMM, MAXPTR, ptr1);
3503 nQMmethod = str_nelem(is->QMmethod, MAXPTR, ptr2);
3504 nQMbasis = str_nelem(is->QMbasis, MAXPTR, ptr3);
3505 if ((nQMmethod != nQMg) || (nQMbasis != nQMg))
3507 gmx_fatal(FARGS, "Invalid QMMM input: %d groups %d basissets"
3508 " and %d methods\n", nQMg, nQMbasis, nQMmethod);
3510 /* group rest, if any, is always MM! */
3511 do_numbering(natoms, groups, nQMg, ptr1, grps, gnames, egcQMMM,
3512 restnm, egrptpALL_GENREST, bVerbose, wi);
3513 nr = nQMg; /*atoms->grps[egcQMMM].nr;*/
3514 ir->opts.ngQM = nQMg;
3515 snew(ir->opts.QMmethod, nr);
3516 snew(ir->opts.QMbasis, nr);
3517 for (i = 0; i < nr; i++)
3519 /* input consists of strings: RHF CASSCF PM3 .. These need to be
3520 * converted to the corresponding enum in names.c
3522 ir->opts.QMmethod[i] = search_QMstring(ptr2[i], eQMmethodNR,
3524 ir->opts.QMbasis[i] = search_QMstring(ptr3[i], eQMbasisNR,
3528 nQMmult = str_nelem(is->QMmult, MAXPTR, ptr1);
3529 nQMcharge = str_nelem(is->QMcharge, MAXPTR, ptr2);
3530 nbSH = str_nelem(is->bSH, MAXPTR, ptr3);
3531 snew(ir->opts.QMmult, nr);
3532 snew(ir->opts.QMcharge, nr);
3533 snew(ir->opts.bSH, nr);
3535 for (i = 0; i < nr; i++)
3537 ir->opts.QMmult[i] = strtol(ptr1[i], NULL, 10);
3538 ir->opts.QMcharge[i] = strtol(ptr2[i], NULL, 10);
3539 ir->opts.bSH[i] = (gmx_strncasecmp(ptr3[i], "Y", 1) == 0);
3542 nCASelec = str_nelem(is->CASelectrons, MAXPTR, ptr1);
3543 nCASorb = str_nelem(is->CASorbitals, MAXPTR, ptr2);
3544 snew(ir->opts.CASelectrons, nr);
3545 snew(ir->opts.CASorbitals, nr);
3546 for (i = 0; i < nr; i++)
3548 ir->opts.CASelectrons[i] = strtol(ptr1[i], NULL, 10);
3549 ir->opts.CASorbitals[i] = strtol(ptr2[i], NULL, 10);
3551 /* special optimization options */
3553 nbOPT = str_nelem(is->bOPT, MAXPTR, ptr1);
3554 nbTS = str_nelem(is->bTS, MAXPTR, ptr2);
3555 snew(ir->opts.bOPT, nr);
3556 snew(ir->opts.bTS, nr);
3557 for (i = 0; i < nr; i++)
3559 ir->opts.bOPT[i] = (gmx_strncasecmp(ptr1[i], "Y", 1) == 0);
3560 ir->opts.bTS[i] = (gmx_strncasecmp(ptr2[i], "Y", 1) == 0);
3562 nSAon = str_nelem(is->SAon, MAXPTR, ptr1);
3563 nSAoff = str_nelem(is->SAoff, MAXPTR, ptr2);
3564 nSAsteps = str_nelem(is->SAsteps, MAXPTR, ptr3);
3565 snew(ir->opts.SAon, nr);
3566 snew(ir->opts.SAoff, nr);
3567 snew(ir->opts.SAsteps, nr);
3569 for (i = 0; i < nr; i++)
3571 ir->opts.SAon[i] = strtod(ptr1[i], NULL);
3572 ir->opts.SAoff[i] = strtod(ptr2[i], NULL);
3573 ir->opts.SAsteps[i] = strtol(ptr3[i], NULL, 10);
3575 /* end of QMMM input */
3579 for (i = 0; (i < egcNR); i++)
3581 fprintf(stderr, "%-16s has %d element(s):", gtypes[i], groups->grps[i].nr);
3582 for (j = 0; (j < groups->grps[i].nr); j++)
3584 fprintf(stderr, " %s", *(groups->grpname[groups->grps[i].nm_ind[j]]));
3586 fprintf(stderr, "\n");
3590 nr = groups->grps[egcENER].nr;
3591 snew(ir->opts.egp_flags, nr*nr);
3593 bExcl = do_egp_flag(ir, groups, "energygrp-excl", is->egpexcl, EGP_EXCL);
3594 if (bExcl && ir->cutoff_scheme == ecutsVERLET)
3596 warning_error(wi, "Energy group exclusions are not (yet) implemented for the Verlet scheme");
3598 if (bExcl && EEL_FULL(ir->coulombtype))
3600 warning(wi, "Can not exclude the lattice Coulomb energy between energy groups");
3603 bTable = do_egp_flag(ir, groups, "energygrp-table", is->egptable, EGP_TABLE);
3604 if (bTable && !(ir->vdwtype == evdwUSER) &&
3605 !(ir->coulombtype == eelUSER) && !(ir->coulombtype == eelPMEUSER) &&
3606 !(ir->coulombtype == eelPMEUSERSWITCH))
3608 gmx_fatal(FARGS, "Can only have energy group pair tables in combination with user tables for VdW and/or Coulomb");
3611 decode_cos(is->efield_x, &(ir->ex[XX]));
3612 decode_cos(is->efield_xt, &(ir->et[XX]));
3613 decode_cos(is->efield_y, &(ir->ex[YY]));
3614 decode_cos(is->efield_yt, &(ir->et[YY]));
3615 decode_cos(is->efield_z, &(ir->ex[ZZ]));
3616 decode_cos(is->efield_zt, &(ir->et[ZZ]));
3620 do_adress_index(ir->adress, groups, gnames, &(ir->opts), wi);
3623 for (i = 0; (i < grps->nr); i++)
3635 static void check_disre(gmx_mtop_t *mtop)
3637 gmx_ffparams_t *ffparams;
3638 t_functype *functype;
3640 int i, ndouble, ftype;
3641 int label, old_label;
3643 if (gmx_mtop_ftype_count(mtop, F_DISRES) > 0)
3645 ffparams = &mtop->ffparams;
3646 functype = ffparams->functype;
3647 ip = ffparams->iparams;
3650 for (i = 0; i < ffparams->ntypes; i++)
3652 ftype = functype[i];
3653 if (ftype == F_DISRES)
3655 label = ip[i].disres.label;
3656 if (label == old_label)
3658 fprintf(stderr, "Distance restraint index %d occurs twice\n", label);
3666 gmx_fatal(FARGS, "Found %d double distance restraint indices,\n"
3667 "probably the parameters for multiple pairs in one restraint "
3668 "are not identical\n", ndouble);
3673 static gmx_bool absolute_reference(t_inputrec *ir, gmx_mtop_t *sys,
3674 gmx_bool posres_only,
3678 gmx_mtop_ilistloop_t iloop;
3688 for (d = 0; d < DIM; d++)
3690 AbsRef[d] = (d < ndof_com(ir) ? 0 : 1);
3692 /* Check for freeze groups */
3693 for (g = 0; g < ir->opts.ngfrz; g++)
3695 for (d = 0; d < DIM; d++)
3697 if (ir->opts.nFreeze[g][d] != 0)
3705 /* Check for position restraints */
3706 iloop = gmx_mtop_ilistloop_init(sys);
3707 while (gmx_mtop_ilistloop_next(iloop, &ilist, &nmol))
3710 (AbsRef[XX] == 0 || AbsRef[YY] == 0 || AbsRef[ZZ] == 0))
3712 for (i = 0; i < ilist[F_POSRES].nr; i += 2)
3714 pr = &sys->ffparams.iparams[ilist[F_POSRES].iatoms[i]];
3715 for (d = 0; d < DIM; d++)
3717 if (pr->posres.fcA[d] != 0)
3723 for (i = 0; i < ilist[F_FBPOSRES].nr; i += 2)
3725 /* Check for flat-bottom posres */
3726 pr = &sys->ffparams.iparams[ilist[F_FBPOSRES].iatoms[i]];
3727 if (pr->fbposres.k != 0)
3729 switch (pr->fbposres.geom)
3731 case efbposresSPHERE:
3732 AbsRef[XX] = AbsRef[YY] = AbsRef[ZZ] = 1;
3734 case efbposresCYLINDER:
3735 AbsRef[XX] = AbsRef[YY] = 1;
3737 case efbposresX: /* d=XX */
3738 case efbposresY: /* d=YY */
3739 case efbposresZ: /* d=ZZ */
3740 d = pr->fbposres.geom - efbposresX;
3744 gmx_fatal(FARGS, " Invalid geometry for flat-bottom position restraint.\n"
3745 "Expected nr between 1 and %d. Found %d\n", efbposresNR-1,
3753 return (AbsRef[XX] != 0 && AbsRef[YY] != 0 && AbsRef[ZZ] != 0);
3757 check_combination_rule_differences(const gmx_mtop_t *mtop, int state,
3758 gmx_bool *bC6ParametersWorkWithGeometricRules,
3759 gmx_bool *bC6ParametersWorkWithLBRules,
3760 gmx_bool *bLBRulesPossible)
3762 int ntypes, tpi, tpj, thisLBdiff, thisgeomdiff;
3765 double geometricdiff, LBdiff;
3766 double c6i, c6j, c12i, c12j;
3767 double c6, c6_geometric, c6_LB;
3768 double sigmai, sigmaj, epsi, epsj;
3769 gmx_bool bCanDoLBRules, bCanDoGeometricRules;
3772 /* A tolerance of 1e-5 seems reasonable for (possibly hand-typed)
3773 * force-field floating point parameters.
3776 ptr = getenv("GMX_LJCOMB_TOL");
3781 sscanf(ptr, "%lf", &dbl);
3785 *bC6ParametersWorkWithLBRules = TRUE;
3786 *bC6ParametersWorkWithGeometricRules = TRUE;
3787 bCanDoLBRules = TRUE;
3788 bCanDoGeometricRules = TRUE;
3789 ntypes = mtop->ffparams.atnr;
3790 snew(typecount, ntypes);
3791 gmx_mtop_count_atomtypes(mtop, state, typecount);
3792 geometricdiff = LBdiff = 0.0;
3793 *bLBRulesPossible = TRUE;
3794 for (tpi = 0; tpi < ntypes; ++tpi)
3796 c6i = mtop->ffparams.iparams[(ntypes + 1) * tpi].lj.c6;
3797 c12i = mtop->ffparams.iparams[(ntypes + 1) * tpi].lj.c12;
3798 for (tpj = tpi; tpj < ntypes; ++tpj)
3800 c6j = mtop->ffparams.iparams[(ntypes + 1) * tpj].lj.c6;
3801 c12j = mtop->ffparams.iparams[(ntypes + 1) * tpj].lj.c12;
3802 c6 = mtop->ffparams.iparams[ntypes * tpi + tpj].lj.c6;
3803 c6_geometric = sqrt(c6i * c6j);
3804 if (!gmx_numzero(c6_geometric))
3806 if (!gmx_numzero(c12i) && !gmx_numzero(c12j))
3808 sigmai = pow(c12i / c6i, 1.0/6.0);
3809 sigmaj = pow(c12j / c6j, 1.0/6.0);
3810 epsi = c6i * c6i /(4.0 * c12i);
3811 epsj = c6j * c6j /(4.0 * c12j);
3812 c6_LB = 4.0 * pow(epsi * epsj, 1.0/2.0) * pow(0.5 * (sigmai + sigmaj), 6);
3816 *bLBRulesPossible = FALSE;
3817 c6_LB = c6_geometric;
3819 bCanDoLBRules = gmx_within_tol(c6_LB, c6, tol);
3822 if (FALSE == bCanDoLBRules)
3824 *bC6ParametersWorkWithLBRules = FALSE;
3827 bCanDoGeometricRules = gmx_within_tol(c6_geometric, c6, tol);
3829 if (FALSE == bCanDoGeometricRules)
3831 *bC6ParametersWorkWithGeometricRules = FALSE;
3839 check_combination_rules(const t_inputrec *ir, const gmx_mtop_t *mtop,
3843 gmx_bool bLBRulesPossible, bC6ParametersWorkWithGeometricRules, bC6ParametersWorkWithLBRules;
3845 check_combination_rule_differences(mtop, 0,
3846 &bC6ParametersWorkWithGeometricRules,
3847 &bC6ParametersWorkWithLBRules,
3849 if (ir->ljpme_combination_rule == eljpmeLB)
3851 if (FALSE == bC6ParametersWorkWithLBRules || FALSE == bLBRulesPossible)
3853 warning(wi, "You are using arithmetic-geometric combination rules "
3854 "in LJ-PME, but your non-bonded C6 parameters do not "
3855 "follow these rules.");
3860 if (FALSE == bC6ParametersWorkWithGeometricRules)
3862 if (ir->eDispCorr != edispcNO)
3864 warning_note(wi, "You are using geometric combination rules in "
3865 "LJ-PME, but your non-bonded C6 parameters do "
3866 "not follow these rules. "
3867 "This will introduce very small errors in the forces and energies in "
3868 "your simulations. Dispersion correction will correct total energy "
3869 "and/or pressure for isotropic systems, but not forces or surface tensions.");
3873 warning_note(wi, "You are using geometric combination rules in "
3874 "LJ-PME, but your non-bonded C6 parameters do "
3875 "not follow these rules. "
3876 "This will introduce very small errors in the forces and energies in "
3877 "your simulations. If your system is homogeneous, consider using dispersion correction "
3878 "for the total energy and pressure.");
3884 void triple_check(const char *mdparin, t_inputrec *ir, gmx_mtop_t *sys,
3888 int i, m, c, nmol, npct;
3889 gmx_bool bCharge, bAcc;
3890 real gdt_max, *mgrp, mt;
3892 gmx_mtop_atomloop_block_t aloopb;
3893 gmx_mtop_atomloop_all_t aloop;
3896 char warn_buf[STRLEN];
3898 set_warning_line(wi, mdparin, -1);
3900 if (EI_DYNAMICS(ir->eI) && !EI_SD(ir->eI) && ir->eI != eiBD &&
3901 ir->comm_mode == ecmNO &&
3902 !(absolute_reference(ir, sys, FALSE, AbsRef) || ir->nsteps <= 10) &&
3903 !ETC_ANDERSEN(ir->etc))
3905 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");
3908 /* Check for pressure coupling with absolute position restraints */
3909 if (ir->epc != epcNO && ir->refcoord_scaling == erscNO)
3911 absolute_reference(ir, sys, TRUE, AbsRef);
3913 for (m = 0; m < DIM; m++)
3915 if (AbsRef[m] && norm2(ir->compress[m]) > 0)
3917 warning(wi, "You are using pressure coupling with absolute position restraints, this will give artifacts. Use the refcoord_scaling option.");
3925 aloopb = gmx_mtop_atomloop_block_init(sys);
3926 while (gmx_mtop_atomloop_block_next(aloopb, &atom, &nmol))
3928 if (atom->q != 0 || atom->qB != 0)
3936 if (EEL_FULL(ir->coulombtype))
3939 "You are using full electrostatics treatment %s for a system without charges.\n"
3940 "This costs a lot of performance for just processing zeros, consider using %s instead.\n",
3941 EELTYPE(ir->coulombtype), EELTYPE(eelCUT));
3942 warning(wi, err_buf);
3947 if (ir->coulombtype == eelCUT && ir->rcoulomb > 0 && !ir->implicit_solvent)
3950 "You are using a plain Coulomb cut-off, which might produce artifacts.\n"
3951 "You might want to consider using %s electrostatics.\n",
3953 warning_note(wi, err_buf);
3957 /* Check if combination rules used in LJ-PME are the same as in the force field */
3958 if (EVDW_PME(ir->vdwtype))
3960 check_combination_rules(ir, sys, wi);
3963 /* Generalized reaction field */
3964 if (ir->opts.ngtc == 0)
3966 sprintf(err_buf, "No temperature coupling while using coulombtype %s",
3968 CHECK(ir->coulombtype == eelGRF);
3972 sprintf(err_buf, "When using coulombtype = %s"
3973 " ref-t for temperature coupling should be > 0",
3975 CHECK((ir->coulombtype == eelGRF) && (ir->opts.ref_t[0] <= 0));
3978 if (ir->eI == eiSD1 &&
3979 (gmx_mtop_ftype_count(sys, F_CONSTR) > 0 ||
3980 gmx_mtop_ftype_count(sys, F_SETTLE) > 0))
3982 sprintf(warn_buf, "With constraints integrator %s is less accurate, consider using %s instead", ei_names[ir->eI], ei_names[eiSD2]);
3983 warning_note(wi, warn_buf);
3987 for (i = 0; (i < sys->groups.grps[egcACC].nr); i++)
3989 for (m = 0; (m < DIM); m++)
3991 if (fabs(ir->opts.acc[i][m]) > 1e-6)
4000 snew(mgrp, sys->groups.grps[egcACC].nr);
4001 aloop = gmx_mtop_atomloop_all_init(sys);
4002 while (gmx_mtop_atomloop_all_next(aloop, &i, &atom))
4004 mgrp[ggrpnr(&sys->groups, egcACC, i)] += atom->m;
4007 for (i = 0; (i < sys->groups.grps[egcACC].nr); i++)
4009 for (m = 0; (m < DIM); m++)
4011 acc[m] += ir->opts.acc[i][m]*mgrp[i];
4015 for (m = 0; (m < DIM); m++)
4017 if (fabs(acc[m]) > 1e-6)
4019 const char *dim[DIM] = { "X", "Y", "Z" };
4021 "Net Acceleration in %s direction, will %s be corrected\n",
4022 dim[m], ir->nstcomm != 0 ? "" : "not");
4023 if (ir->nstcomm != 0 && m < ndof_com(ir))
4026 for (i = 0; (i < sys->groups.grps[egcACC].nr); i++)
4028 ir->opts.acc[i][m] -= acc[m];
4036 if (ir->efep != efepNO && ir->fepvals->sc_alpha != 0 &&
4037 !gmx_within_tol(sys->ffparams.reppow, 12.0, 10*GMX_DOUBLE_EPS))
4039 gmx_fatal(FARGS, "Soft-core interactions are only supported with VdW repulsion power 12");
4042 if (ir->ePull != epullNO)
4044 gmx_bool bPullAbsoluteRef;
4046 bPullAbsoluteRef = FALSE;
4047 for (i = 0; i < ir->pull->ncoord; i++)
4049 bPullAbsoluteRef = bPullAbsoluteRef ||
4050 ir->pull->coord[i].group[0] == 0 ||
4051 ir->pull->coord[i].group[1] == 0;
4053 if (bPullAbsoluteRef)
4055 absolute_reference(ir, sys, FALSE, AbsRef);
4056 for (m = 0; m < DIM; m++)
4058 if (ir->pull->dim[m] && !AbsRef[m])
4060 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.");
4066 if (ir->pull->eGeom == epullgDIRPBC)
4068 for (i = 0; i < 3; i++)
4070 for (m = 0; m <= i; m++)
4072 if ((ir->epc != epcNO && ir->compress[i][m] != 0) ||
4073 ir->deform[i][m] != 0)
4075 for (c = 0; c < ir->pull->ncoord; c++)
4077 if (ir->pull->coord[c].vec[m] != 0)
4079 gmx_fatal(FARGS, "Can not have dynamic box while using pull geometry '%s' (dim %c)", EPULLGEOM(ir->pull->eGeom), 'x'+m);
4091 void double_check(t_inputrec *ir, matrix box, gmx_bool bConstr, warninp_t wi)
4095 char warn_buf[STRLEN];
4098 ptr = check_box(ir->ePBC, box);
4101 warning_error(wi, ptr);
4104 if (bConstr && ir->eConstrAlg == econtSHAKE)
4106 if (ir->shake_tol <= 0.0)
4108 sprintf(warn_buf, "ERROR: shake-tol must be > 0 instead of %g\n",
4110 warning_error(wi, warn_buf);
4113 if (IR_TWINRANGE(*ir) && ir->nstlist > 1)
4115 sprintf(warn_buf, "With twin-range cut-off's and SHAKE the virial and the pressure are incorrect.");
4116 if (ir->epc == epcNO)
4118 warning(wi, warn_buf);
4122 warning_error(wi, warn_buf);
4127 if ( (ir->eConstrAlg == econtLINCS) && bConstr)
4129 /* If we have Lincs constraints: */
4130 if (ir->eI == eiMD && ir->etc == etcNO &&
4131 ir->eConstrAlg == econtLINCS && ir->nLincsIter == 1)
4133 sprintf(warn_buf, "For energy conservation with LINCS, lincs_iter should be 2 or larger.\n");
4134 warning_note(wi, warn_buf);
4137 if ((ir->eI == eiCG || ir->eI == eiLBFGS) && (ir->nProjOrder < 8))
4139 sprintf(warn_buf, "For accurate %s with LINCS constraints, lincs-order should be 8 or more.", ei_names[ir->eI]);
4140 warning_note(wi, warn_buf);
4142 if (ir->epc == epcMTTK)
4144 warning_error(wi, "MTTK not compatible with lincs -- use shake instead.");
4148 if (ir->LincsWarnAngle > 90.0)
4150 sprintf(warn_buf, "lincs-warnangle can not be larger than 90 degrees, setting it to 90.\n");
4151 warning(wi, warn_buf);
4152 ir->LincsWarnAngle = 90.0;
4155 if (ir->ePBC != epbcNONE)
4157 if (ir->nstlist == 0)
4159 warning(wi, "With nstlist=0 atoms are only put into the box at step 0, therefore drifting atoms might cause the simulation to crash.");
4161 bTWIN = (ir->rlistlong > ir->rlist);
4162 if (ir->ns_type == ensGRID)
4164 if (sqr(ir->rlistlong) >= max_cutoff2(ir->ePBC, box))
4166 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",
4167 bTWIN ? (ir->rcoulomb == ir->rlistlong ? "rcoulomb" : "rvdw") : "rlist");
4168 warning_error(wi, warn_buf);
4173 min_size = min(box[XX][XX], min(box[YY][YY], box[ZZ][ZZ]));
4174 if (2*ir->rlistlong >= min_size)
4176 sprintf(warn_buf, "ERROR: One of the box lengths is smaller than twice the cut-off length. Increase the box size or decrease rlist.");
4177 warning_error(wi, warn_buf);
4180 fprintf(stderr, "Grid search might allow larger cut-off's than simple search with triclinic boxes.");
4187 void check_chargegroup_radii(const gmx_mtop_t *mtop, const t_inputrec *ir,
4191 real rvdw1, rvdw2, rcoul1, rcoul2;
4192 char warn_buf[STRLEN];
4194 calc_chargegroup_radii(mtop, x, &rvdw1, &rvdw2, &rcoul1, &rcoul2);
4198 printf("Largest charge group radii for Van der Waals: %5.3f, %5.3f nm\n",
4203 printf("Largest charge group radii for Coulomb: %5.3f, %5.3f nm\n",
4209 if (rvdw1 + rvdw2 > ir->rlist ||
4210 rcoul1 + rcoul2 > ir->rlist)
4213 "The sum of the two largest charge group radii (%f) "
4214 "is larger than rlist (%f)\n",
4215 max(rvdw1+rvdw2, rcoul1+rcoul2), ir->rlist);
4216 warning(wi, warn_buf);
4220 /* Here we do not use the zero at cut-off macro,
4221 * since user defined interactions might purposely
4222 * not be zero at the cut-off.
4224 if (ir_vdw_is_zero_at_cutoff(ir) &&
4225 rvdw1 + rvdw2 > ir->rlistlong - ir->rvdw)
4227 sprintf(warn_buf, "The sum of the two largest charge group "
4228 "radii (%f) is larger than %s (%f) - rvdw (%f).\n"
4229 "With exact cut-offs, better performance can be "
4230 "obtained with cutoff-scheme = %s, because it "
4231 "does not use charge groups at all.",
4233 ir->rlistlong > ir->rlist ? "rlistlong" : "rlist",
4234 ir->rlistlong, ir->rvdw,
4235 ecutscheme_names[ecutsVERLET]);
4238 warning(wi, warn_buf);
4242 warning_note(wi, warn_buf);
4245 if (ir_coulomb_is_zero_at_cutoff(ir) &&
4246 rcoul1 + rcoul2 > ir->rlistlong - ir->rcoulomb)
4248 sprintf(warn_buf, "The sum of the two largest charge group radii (%f) is larger than %s (%f) - rcoulomb (%f).\n"
4249 "With exact cut-offs, better performance can be obtained with cutoff-scheme = %s, because it does not use charge groups at all.",
4251 ir->rlistlong > ir->rlist ? "rlistlong" : "rlist",
4252 ir->rlistlong, ir->rcoulomb,
4253 ecutscheme_names[ecutsVERLET]);
4256 warning(wi, warn_buf);
4260 warning_note(wi, warn_buf);