1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
4 * This source code is part of
8 * GROningen MAchine for Chemical Simulations
11 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
12 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
13 * Copyright (c) 2001-2004, The GROMACS development team,
14 * check out http://www.gromacs.org for more information.
16 * This program is free software; you can redistribute it and/or
17 * modify it under the terms of the GNU General Public License
18 * as published by the Free Software Foundation; either version 2
19 * of the License, or (at your option) any later version.
21 * If you want to redistribute modifications, please consider that
22 * scientific software is very special. Version control is crucial -
23 * bugs must be traceable. We will be happy to consider code for
24 * inclusion in the official distribution, but derived work must not
25 * be called official GROMACS. Details are found in the README & COPYING
26 * files - if they are missing, get the official version at www.gromacs.org.
28 * To help us fund GROMACS development, we humbly ask that you cite
29 * the papers on the package - you can find them in the top README file.
31 * For more info, check our website at http://www.gromacs.org
34 * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
48 #include "gmx_fatal.h"
61 #include "mtop_util.h"
62 #include "chargegroup.h"
67 #define MAXLAMBDAS 1024
69 /* Resource parameters
70 * Do not change any of these until you read the instruction
71 * in readinp.h. Some cpp's do not take spaces after the backslash
72 * (like the c-shell), which will give you a very weird compiler
76 static char tcgrps[STRLEN], tau_t[STRLEN], ref_t[STRLEN],
77 acc[STRLEN], accgrps[STRLEN], freeze[STRLEN], frdim[STRLEN],
78 energy[STRLEN], user1[STRLEN], user2[STRLEN], vcm[STRLEN], xtc_grps[STRLEN],
79 couple_moltype[STRLEN], orirefitgrp[STRLEN], egptable[STRLEN], egpexcl[STRLEN],
80 wall_atomtype[STRLEN], wall_density[STRLEN], deform[STRLEN], QMMM[STRLEN];
81 static char fep_lambda[efptNR][STRLEN];
82 static char lambda_weights[STRLEN];
83 static char **pull_grp;
84 static char **rot_grp;
85 static char anneal[STRLEN], anneal_npoints[STRLEN],
86 anneal_time[STRLEN], anneal_temp[STRLEN];
87 static char QMmethod[STRLEN], QMbasis[STRLEN], QMcharge[STRLEN], QMmult[STRLEN],
88 bSH[STRLEN], CASorbitals[STRLEN], CASelectrons[STRLEN], SAon[STRLEN],
89 SAoff[STRLEN], SAsteps[STRLEN], bTS[STRLEN], bOPT[STRLEN];
90 static char efield_x[STRLEN], efield_xt[STRLEN], efield_y[STRLEN],
91 efield_yt[STRLEN], efield_z[STRLEN], efield_zt[STRLEN];
94 egrptpALL, /* All particles have to be a member of a group. */
95 egrptpALL_GENREST, /* A rest group with name is generated for particles *
96 * that are not part of any group. */
97 egrptpPART, /* As egrptpALL_GENREST, but no name is generated *
98 * for the rest group. */
99 egrptpONE /* Merge all selected groups into one group, *
100 * make a rest group for the remaining particles. */
103 static const char *constraints[eshNR+1] = {
104 "none", "h-bonds", "all-bonds", "h-angles", "all-angles", NULL
107 static const char *couple_lam[ecouplamNR+1] = {
108 "vdw-q", "vdw", "q", "none", NULL
111 void init_ir(t_inputrec *ir, t_gromppopts *opts)
113 snew(opts->include, STRLEN);
114 snew(opts->define, STRLEN);
115 snew(ir->fepvals, 1);
116 snew(ir->expandedvals, 1);
117 snew(ir->simtempvals, 1);
120 static void GetSimTemps(int ntemps, t_simtemp *simtemp, double *temperature_lambdas)
125 for (i = 0; i < ntemps; i++)
127 /* simple linear scaling -- allows more control */
128 if (simtemp->eSimTempScale == esimtempLINEAR)
130 simtemp->temperatures[i] = simtemp->simtemp_low + (simtemp->simtemp_high-simtemp->simtemp_low)*temperature_lambdas[i];
132 else if (simtemp->eSimTempScale == esimtempGEOMETRIC) /* should give roughly equal acceptance for constant heat capacity . . . */
134 simtemp->temperatures[i] = simtemp->simtemp_low * pow(simtemp->simtemp_high/simtemp->simtemp_low, (1.0*i)/(ntemps-1));
136 else if (simtemp->eSimTempScale == esimtempEXPONENTIAL)
138 simtemp->temperatures[i] = simtemp->simtemp_low + (simtemp->simtemp_high-simtemp->simtemp_low)*((exp(temperature_lambdas[i])-1)/(exp(1.0)-1));
143 sprintf(errorstr, "eSimTempScale=%d not defined", simtemp->eSimTempScale);
144 gmx_fatal(FARGS, errorstr);
151 static void _low_check(gmx_bool b, char *s, warninp_t wi)
155 warning_error(wi, s);
159 static void check_nst(const char *desc_nst, int nst,
160 const char *desc_p, int *p,
165 if (*p > 0 && *p % nst != 0)
167 /* Round up to the next multiple of nst */
168 *p = ((*p)/nst + 1)*nst;
169 sprintf(buf, "%s should be a multiple of %s, changing %s to %d\n",
170 desc_p, desc_nst, desc_p, *p);
175 static gmx_bool ir_NVE(const t_inputrec *ir)
177 return ((ir->eI == eiMD || EI_VV(ir->eI)) && ir->etc == etcNO);
180 static int lcd(int n1, int n2)
185 for (i = 2; (i <= n1 && i <= n2); i++)
187 if (n1 % i == 0 && n2 % i == 0)
196 static void process_interaction_modifier(const t_inputrec *ir, int *eintmod)
198 if (*eintmod == eintmodPOTSHIFT_VERLET)
200 if (ir->cutoff_scheme == ecutsVERLET)
202 *eintmod = eintmodPOTSHIFT;
206 *eintmod = eintmodNONE;
211 void check_ir(const char *mdparin, t_inputrec *ir, t_gromppopts *opts,
213 /* Check internal consistency */
215 /* Strange macro: first one fills the err_buf, and then one can check
216 * the condition, which will print the message and increase the error
219 #define CHECK(b) _low_check(b, err_buf, wi)
220 char err_buf[256], warn_buf[STRLEN];
226 t_lambda *fep = ir->fepvals;
227 t_expanded *expand = ir->expandedvals;
229 set_warning_line(wi, mdparin, -1);
231 /* BASIC CUT-OFF STUFF */
232 if (ir->rcoulomb < 0)
234 warning_error(wi, "rcoulomb should be >= 0");
238 warning_error(wi, "rvdw should be >= 0");
241 !(ir->cutoff_scheme == ecutsVERLET && ir->verletbuf_drift > 0))
243 warning_error(wi, "rlist should be >= 0");
246 process_interaction_modifier(ir, &ir->coulomb_modifier);
247 process_interaction_modifier(ir, &ir->vdw_modifier);
249 if (ir->cutoff_scheme == ecutsGROUP)
251 /* BASIC CUT-OFF STUFF */
252 if (ir->rlist == 0 ||
253 !((EEL_MIGHT_BE_ZERO_AT_CUTOFF(ir->coulombtype) && ir->rcoulomb > ir->rlist) ||
254 (EVDW_MIGHT_BE_ZERO_AT_CUTOFF(ir->vdwtype) && ir->rvdw > ir->rlist)))
256 /* No switched potential and/or no twin-range:
257 * we can set the long-range cut-off to the maximum of the other cut-offs.
259 ir->rlistlong = max_cutoff(ir->rlist, max_cutoff(ir->rvdw, ir->rcoulomb));
261 else if (ir->rlistlong < 0)
263 ir->rlistlong = max_cutoff(ir->rlist, max_cutoff(ir->rvdw, ir->rcoulomb));
264 sprintf(warn_buf, "rlistlong was not set, setting it to %g (no buffer)",
266 warning(wi, warn_buf);
268 if (ir->rlistlong == 0 && ir->ePBC != epbcNONE)
270 warning_error(wi, "Can not have an infinite cut-off with PBC");
272 if (ir->rlistlong > 0 && (ir->rlist == 0 || ir->rlistlong < ir->rlist))
274 warning_error(wi, "rlistlong can not be shorter than rlist");
276 if (IR_TWINRANGE(*ir) && ir->nstlist <= 0)
278 warning_error(wi, "Can not have nstlist<=0 with twin-range interactions");
282 if (ir->rlistlong == ir->rlist)
286 else if (ir->rlistlong > ir->rlist && ir->nstcalclr == 0)
288 warning_error(wi, "With different cutoffs for electrostatics and VdW, nstcalclr must be -1 or a positive number");
291 if (ir->cutoff_scheme == ecutsVERLET)
295 /* Normal Verlet type neighbor-list, currently only limited feature support */
296 if (inputrec2nboundeddim(ir) < 3)
298 warning_error(wi, "With Verlet lists only full pbc or pbc=xy with walls is supported");
300 if (ir->rcoulomb != ir->rvdw)
302 warning_error(wi, "With Verlet lists rcoulomb!=rvdw is not supported");
304 if (ir->vdwtype != evdwCUT)
306 warning_error(wi, "With Verlet lists only cut-off LJ interactions are supported");
308 if (!(ir->coulombtype == eelCUT ||
309 (EEL_RF(ir->coulombtype) && ir->coulombtype != eelRF_NEC) ||
310 EEL_PME(ir->coulombtype) || ir->coulombtype == eelEWALD))
312 warning_error(wi, "With Verlet lists only cut-off, reaction-field, PME and Ewald electrostatics are supported");
315 if (ir->nstlist <= 0)
317 warning_error(wi, "With Verlet lists nstlist should be larger than 0");
320 if (ir->nstlist < 10)
322 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.");
325 rc_max = max(ir->rvdw, ir->rcoulomb);
327 if (ir->verletbuf_drift <= 0)
329 if (ir->verletbuf_drift == 0)
331 warning_error(wi, "Can not have an energy drift of exactly 0");
334 if (ir->rlist < rc_max)
336 warning_error(wi, "With verlet lists rlist can not be smaller than rvdw or rcoulomb");
339 if (ir->rlist == rc_max && ir->nstlist > 1)
341 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.");
346 if (ir->rlist > rc_max)
348 warning_note(wi, "You have set rlist larger than the interaction cut-off, but you also have verlet-buffer-drift > 0. Will set rlist using verlet-buffer-drift.");
351 if (ir->nstlist == 1)
353 /* No buffer required */
358 if (EI_DYNAMICS(ir->eI))
360 if (EI_MD(ir->eI) && ir->etc == etcNO)
362 warning_error(wi, "Temperature coupling is required for calculating rlist using the energy drift with verlet-buffer-drift > 0. Either use temperature coupling or set rlist yourself together with verlet-buffer-drift = -1.");
365 if (inputrec2nboundeddim(ir) < 3)
367 warning_error(wi, "The box volume is required for calculating rlist from the energy drift with verlet-buffer-drift > 0. You are using at least one unbounded dimension, so no volume can be computed. Either use a finite box, or set rlist yourself together with verlet-buffer-drift = -1.");
369 /* Set rlist temporarily so we can continue processing */
374 /* Set the buffer to 5% of the cut-off */
375 ir->rlist = 1.05*rc_max;
380 /* No twin-range calculations with Verlet lists */
381 ir->rlistlong = ir->rlist;
384 if (ir->nstcalclr == -1)
386 /* if rlist=rlistlong, this will later be changed to nstcalclr=0 */
387 ir->nstcalclr = ir->nstlist;
389 else if (ir->nstcalclr > 0)
391 if (ir->nstlist > 0 && (ir->nstlist % ir->nstcalclr != 0))
393 warning_error(wi, "nstlist must be evenly divisible by nstcalclr. Use nstcalclr = -1 to automatically follow nstlist");
396 else if (ir->nstcalclr < -1)
398 warning_error(wi, "nstcalclr must be a positive number (divisor of nstcalclr), or -1 to follow nstlist.");
401 if (EEL_PME(ir->coulombtype) && ir->rcoulomb > ir->rvdw && ir->nstcalclr > 1)
403 warning_error(wi, "When used with PME, the long-range component of twin-range interactions must be updated every step (nstcalclr)");
406 /* GENERAL INTEGRATOR STUFF */
407 if (!(ir->eI == eiMD || EI_VV(ir->eI)))
411 if (ir->eI == eiVVAK)
413 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]);
414 warning_note(wi, warn_buf);
416 if (!EI_DYNAMICS(ir->eI))
420 if (EI_DYNAMICS(ir->eI))
422 if (ir->nstcalcenergy < 0)
424 ir->nstcalcenergy = ir_optimal_nstcalcenergy(ir);
425 if (ir->nstenergy != 0 && ir->nstenergy < ir->nstcalcenergy)
427 /* nstcalcenergy larger than nstener does not make sense.
428 * We ideally want nstcalcenergy=nstener.
432 ir->nstcalcenergy = lcd(ir->nstenergy, ir->nstlist);
436 ir->nstcalcenergy = ir->nstenergy;
440 else if ( (ir->nstenergy > 0 && ir->nstcalcenergy > ir->nstenergy) ||
441 (ir->efep != efepNO && ir->fepvals->nstdhdl > 0 &&
442 (ir->nstcalcenergy > ir->fepvals->nstdhdl) ) )
445 const char *nsten = "nstenergy";
446 const char *nstdh = "nstdhdl";
447 const char *min_name = nsten;
448 int min_nst = ir->nstenergy;
450 /* find the smallest of ( nstenergy, nstdhdl ) */
451 if (ir->efep != efepNO && ir->fepvals->nstdhdl > 0 &&
452 (ir->fepvals->nstdhdl < ir->nstenergy) )
454 min_nst = ir->fepvals->nstdhdl;
457 /* If the user sets nstenergy small, we should respect that */
459 "Setting nstcalcenergy (%d) equal to %s (%d)",
460 ir->nstcalcenergy, min_name, min_nst);
461 warning_note(wi, warn_buf);
462 ir->nstcalcenergy = min_nst;
465 if (ir->epc != epcNO)
467 if (ir->nstpcouple < 0)
469 ir->nstpcouple = ir_optimal_nstpcouple(ir);
472 if (IR_TWINRANGE(*ir))
474 check_nst("nstlist", ir->nstlist,
475 "nstcalcenergy", &ir->nstcalcenergy, wi);
476 if (ir->epc != epcNO)
478 check_nst("nstlist", ir->nstlist,
479 "nstpcouple", &ir->nstpcouple, wi);
483 if (ir->nstcalcenergy > 0)
485 if (ir->efep != efepNO)
487 /* nstdhdl should be a multiple of nstcalcenergy */
488 check_nst("nstcalcenergy", ir->nstcalcenergy,
489 "nstdhdl", &ir->fepvals->nstdhdl, wi);
490 /* nstexpanded should be a multiple of nstcalcenergy */
491 check_nst("nstcalcenergy", ir->nstcalcenergy,
492 "nstexpanded", &ir->expandedvals->nstexpanded, wi);
494 /* for storing exact averages nstenergy should be
495 * a multiple of nstcalcenergy
497 check_nst("nstcalcenergy", ir->nstcalcenergy,
498 "nstenergy", &ir->nstenergy, wi);
503 if ((EI_SD(ir->eI) || ir->eI == eiBD) &&
504 ir->bContinuation && ir->ld_seed != -1)
506 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)");
512 sprintf(err_buf, "TPI only works with pbc = %s", epbc_names[epbcXYZ]);
513 CHECK(ir->ePBC != epbcXYZ);
514 sprintf(err_buf, "TPI only works with ns = %s", ens_names[ensGRID]);
515 CHECK(ir->ns_type != ensGRID);
516 sprintf(err_buf, "with TPI nstlist should be larger than zero");
517 CHECK(ir->nstlist <= 0);
518 sprintf(err_buf, "TPI does not work with full electrostatics other than PME");
519 CHECK(EEL_FULL(ir->coulombtype) && !EEL_PME(ir->coulombtype));
523 if ( (opts->nshake > 0) && (opts->bMorse) )
526 "Using morse bond-potentials while constraining bonds is useless");
527 warning(wi, warn_buf);
530 if ((EI_SD(ir->eI) || ir->eI == eiBD) &&
531 ir->bContinuation && ir->ld_seed != -1)
533 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)");
535 /* verify simulated tempering options */
539 gmx_bool bAllTempZero = TRUE;
540 for (i = 0; i < fep->n_lambda; i++)
542 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]);
543 CHECK((fep->all_lambda[efptTEMPERATURE][i] < 0) || (fep->all_lambda[efptTEMPERATURE][i] > 1));
544 if (fep->all_lambda[efptTEMPERATURE][i] > 0)
546 bAllTempZero = FALSE;
549 sprintf(err_buf, "if simulated tempering is on, temperature-lambdas may not be all zero");
550 CHECK(bAllTempZero == TRUE);
552 sprintf(err_buf, "Simulated tempering is currently only compatible with md-vv");
553 CHECK(ir->eI != eiVV);
555 /* check compatability of the temperature coupling with simulated tempering */
557 if (ir->etc == etcNOSEHOOVER)
559 sprintf(warn_buf, "Nose-Hoover based temperature control such as [%s] my not be entirelyconsistent with simulated tempering", etcoupl_names[ir->etc]);
560 warning_note(wi, warn_buf);
563 /* check that the temperatures make sense */
565 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);
566 CHECK(ir->simtempvals->simtemp_high <= ir->simtempvals->simtemp_low);
568 sprintf(err_buf, "Higher simulated tempering temperature (%g) must be >= zero", ir->simtempvals->simtemp_high);
569 CHECK(ir->simtempvals->simtemp_high <= 0);
571 sprintf(err_buf, "Lower simulated tempering temperature (%g) must be >= zero", ir->simtempvals->simtemp_low);
572 CHECK(ir->simtempvals->simtemp_low <= 0);
575 /* verify free energy options */
577 if (ir->efep != efepNO)
580 sprintf(err_buf, "The soft-core power is %d and can only be 1 or 2",
582 CHECK(fep->sc_alpha != 0 && fep->sc_power != 1 && fep->sc_power != 2);
584 sprintf(err_buf, "The soft-core sc-r-power is %d and can only be 6 or 48",
585 (int)fep->sc_r_power);
586 CHECK(fep->sc_alpha != 0 && fep->sc_r_power != 6.0 && fep->sc_r_power != 48.0);
588 sprintf(err_buf, "Can't use postive delta-lambda (%g) if initial state/lambda does not start at zero", fep->delta_lambda);
589 CHECK(fep->delta_lambda > 0 && ((fep->init_fep_state > 0) || (fep->init_lambda > 0)));
591 sprintf(err_buf, "Can't use postive delta-lambda (%g) with expanded ensemble simulations", fep->delta_lambda);
592 CHECK(fep->delta_lambda > 0 && (ir->efep == efepEXPANDED));
594 sprintf(err_buf, "Can only use expanded ensemble with md-vv for now; should be supported for other integrators in 5.0");
595 CHECK(!(EI_VV(ir->eI)) && (ir->efep == efepEXPANDED));
597 sprintf(err_buf, "Free-energy not implemented for Ewald");
598 CHECK(ir->coulombtype == eelEWALD);
600 /* check validty of lambda inputs */
601 if (fep->n_lambda == 0)
603 /* Clear output in case of no states:*/
604 sprintf(err_buf, "init-lambda-state set to %d: no lambda states are defined.", fep->init_fep_state);
605 CHECK((fep->init_fep_state >= 0) && (fep->n_lambda == 0));
609 sprintf(err_buf, "initial thermodynamic state %d does not exist, only goes to %d", fep->init_fep_state, fep->n_lambda-1);
610 CHECK((fep->init_fep_state >= fep->n_lambda));
613 sprintf(err_buf, "Lambda state must be set, either with init-lambda-state or with init-lambda");
614 CHECK((fep->init_fep_state < 0) && (fep->init_lambda < 0));
616 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",
617 fep->init_lambda, fep->init_fep_state);
618 CHECK((fep->init_fep_state >= 0) && (fep->init_lambda >= 0));
622 if ((fep->init_lambda >= 0) && (fep->delta_lambda == 0))
626 for (i = 0; i < efptNR; i++)
628 if (fep->separate_dvdl[i])
633 if (n_lambda_terms > 1)
635 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.");
636 warning(wi, warn_buf);
639 if (n_lambda_terms < 2 && fep->n_lambda > 0)
642 "init-lambda is deprecated for setting lambda state (except for slow growth). Use init-lambda-state instead.");
646 for (j = 0; j < efptNR; j++)
648 for (i = 0; i < fep->n_lambda; i++)
650 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]);
651 CHECK((fep->all_lambda[j][i] < 0) || (fep->all_lambda[j][i] > 1));
655 if ((fep->sc_alpha > 0) && (!fep->bScCoul))
657 for (i = 0; i < fep->n_lambda; i++)
659 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],
660 fep->all_lambda[efptCOUL][i]);
661 CHECK((fep->sc_alpha > 0) &&
662 (((fep->all_lambda[efptCOUL][i] > 0.0) &&
663 (fep->all_lambda[efptCOUL][i] < 1.0)) &&
664 ((fep->all_lambda[efptVDW][i] > 0.0) &&
665 (fep->all_lambda[efptVDW][i] < 1.0))));
669 if ((fep->bScCoul) && (EEL_PME(ir->coulombtype)))
671 real sigma, lambda, r_sc;
674 /* Maximum estimate for A and B charges equal with lambda power 1 */
676 r_sc = pow(lambda*fep->sc_alpha*pow(sigma/ir->rcoulomb, fep->sc_r_power) + 1.0, 1.0/fep->sc_r_power);
677 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.",
679 sigma, lambda, r_sc - 1.0, ir->ewald_rtol);
680 warning_note(wi, warn_buf);
683 /* Free Energy Checks -- In an ideal world, slow growth and FEP would
684 be treated differently, but that's the next step */
686 for (i = 0; i < efptNR; i++)
688 for (j = 0; j < fep->n_lambda; j++)
690 sprintf(err_buf, "%s[%d] must be between 0 and 1", efpt_names[i], j);
691 CHECK((fep->all_lambda[i][j] < 0) || (fep->all_lambda[i][j] > 1));
696 if ((ir->bSimTemp) || (ir->efep == efepEXPANDED))
699 expand = ir->expandedvals;
701 /* checking equilibration of weights inputs for validity */
703 sprintf(err_buf, "weight-equil-number-all-lambda (%d) is ignored if lmc-weights-equil is not equal to %s",
704 expand->equil_n_at_lam, elmceq_names[elmceqNUMATLAM]);
705 CHECK((expand->equil_n_at_lam > 0) && (expand->elmceq != elmceqNUMATLAM));
707 sprintf(err_buf, "weight-equil-number-samples (%d) is ignored if lmc-weights-equil is not equal to %s",
708 expand->equil_samples, elmceq_names[elmceqSAMPLES]);
709 CHECK((expand->equil_samples > 0) && (expand->elmceq != elmceqSAMPLES));
711 sprintf(err_buf, "weight-equil-number-steps (%d) is ignored if lmc-weights-equil is not equal to %s",
712 expand->equil_steps, elmceq_names[elmceqSTEPS]);
713 CHECK((expand->equil_steps > 0) && (expand->elmceq != elmceqSTEPS));
715 sprintf(err_buf, "weight-equil-wl-delta (%d) is ignored if lmc-weights-equil is not equal to %s",
716 expand->equil_samples, elmceq_names[elmceqWLDELTA]);
717 CHECK((expand->equil_wl_delta > 0) && (expand->elmceq != elmceqWLDELTA));
719 sprintf(err_buf, "weight-equil-count-ratio (%f) is ignored if lmc-weights-equil is not equal to %s",
720 expand->equil_ratio, elmceq_names[elmceqRATIO]);
721 CHECK((expand->equil_ratio > 0) && (expand->elmceq != elmceqRATIO));
723 sprintf(err_buf, "weight-equil-number-all-lambda (%d) must be a positive integer if lmc-weights-equil=%s",
724 expand->equil_n_at_lam, elmceq_names[elmceqNUMATLAM]);
725 CHECK((expand->equil_n_at_lam <= 0) && (expand->elmceq == elmceqNUMATLAM));
727 sprintf(err_buf, "weight-equil-number-samples (%d) must be a positive integer if lmc-weights-equil=%s",
728 expand->equil_samples, elmceq_names[elmceqSAMPLES]);
729 CHECK((expand->equil_samples <= 0) && (expand->elmceq == elmceqSAMPLES));
731 sprintf(err_buf, "weight-equil-number-steps (%d) must be a positive integer if lmc-weights-equil=%s",
732 expand->equil_steps, elmceq_names[elmceqSTEPS]);
733 CHECK((expand->equil_steps <= 0) && (expand->elmceq == elmceqSTEPS));
735 sprintf(err_buf, "weight-equil-wl-delta (%f) must be > 0 if lmc-weights-equil=%s",
736 expand->equil_wl_delta, elmceq_names[elmceqWLDELTA]);
737 CHECK((expand->equil_wl_delta <= 0) && (expand->elmceq == elmceqWLDELTA));
739 sprintf(err_buf, "weight-equil-count-ratio (%f) must be > 0 if lmc-weights-equil=%s",
740 expand->equil_ratio, elmceq_names[elmceqRATIO]);
741 CHECK((expand->equil_ratio <= 0) && (expand->elmceq == elmceqRATIO));
743 sprintf(err_buf, "lmc-weights-equil=%s only possible when lmc-stats = %s or lmc-stats %s",
744 elmceq_names[elmceqWLDELTA], elamstats_names[elamstatsWL], elamstats_names[elamstatsWWL]);
745 CHECK((expand->elmceq == elmceqWLDELTA) && (!EWL(expand->elamstats)));
747 sprintf(err_buf, "lmc-repeats (%d) must be greater than 0", expand->lmc_repeats);
748 CHECK((expand->lmc_repeats <= 0));
749 sprintf(err_buf, "minimum-var-min (%d) must be greater than 0", expand->minvarmin);
750 CHECK((expand->minvarmin <= 0));
751 sprintf(err_buf, "weight-c-range (%d) must be greater or equal to 0", expand->c_range);
752 CHECK((expand->c_range < 0));
753 sprintf(err_buf, "init-lambda-state (%d) must be zero if lmc-forced-nstart (%d)> 0 and lmc-move != 'no'",
754 fep->init_fep_state, expand->lmc_forced_nstart);
755 CHECK((fep->init_fep_state != 0) && (expand->lmc_forced_nstart > 0) && (expand->elmcmove != elmcmoveNO));
756 sprintf(err_buf, "lmc-forced-nstart (%d) must not be negative", expand->lmc_forced_nstart);
757 CHECK((expand->lmc_forced_nstart < 0));
758 sprintf(err_buf, "init-lambda-state (%d) must be in the interval [0,number of lambdas)", fep->init_fep_state);
759 CHECK((fep->init_fep_state < 0) || (fep->init_fep_state >= fep->n_lambda));
761 sprintf(err_buf, "init-wl-delta (%f) must be greater than or equal to 0", expand->init_wl_delta);
762 CHECK((expand->init_wl_delta < 0));
763 sprintf(err_buf, "wl-ratio (%f) must be between 0 and 1", expand->wl_ratio);
764 CHECK((expand->wl_ratio <= 0) || (expand->wl_ratio >= 1));
765 sprintf(err_buf, "wl-scale (%f) must be between 0 and 1", expand->wl_scale);
766 CHECK((expand->wl_scale <= 0) || (expand->wl_scale >= 1));
768 /* if there is no temperature control, we need to specify an MC temperature */
769 sprintf(err_buf, "If there is no temperature control, and lmc-mcmove!= 'no',mc_temperature must be set to a positive number");
770 if (expand->nstTij > 0)
772 sprintf(err_buf, "nst-transition-matrix (%d) must be an integer multiple of nstlog (%d)",
773 expand->nstTij, ir->nstlog);
774 CHECK((mod(expand->nstTij, ir->nstlog) != 0));
779 sprintf(err_buf, "walls only work with pbc=%s", epbc_names[epbcXY]);
780 CHECK(ir->nwall && ir->ePBC != epbcXY);
783 if (ir->ePBC != epbcXYZ && ir->nwall != 2)
785 if (ir->ePBC == epbcNONE)
787 if (ir->epc != epcNO)
789 warning(wi, "Turning off pressure coupling for vacuum system");
795 sprintf(err_buf, "Can not have pressure coupling with pbc=%s",
796 epbc_names[ir->ePBC]);
797 CHECK(ir->epc != epcNO);
799 sprintf(err_buf, "Can not have Ewald with pbc=%s", epbc_names[ir->ePBC]);
800 CHECK(EEL_FULL(ir->coulombtype));
802 sprintf(err_buf, "Can not have dispersion correction with pbc=%s",
803 epbc_names[ir->ePBC]);
804 CHECK(ir->eDispCorr != edispcNO);
807 if (ir->rlist == 0.0)
809 sprintf(err_buf, "can only have neighborlist cut-off zero (=infinite)\n"
810 "with coulombtype = %s or coulombtype = %s\n"
811 "without periodic boundary conditions (pbc = %s) and\n"
812 "rcoulomb and rvdw set to zero",
813 eel_names[eelCUT], eel_names[eelUSER], epbc_names[epbcNONE]);
814 CHECK(((ir->coulombtype != eelCUT) && (ir->coulombtype != eelUSER)) ||
815 (ir->ePBC != epbcNONE) ||
816 (ir->rcoulomb != 0.0) || (ir->rvdw != 0.0));
820 warning_error(wi, "Can not have heuristic neighborlist updates without cut-off");
824 warning_note(wi, "Simulating without cut-offs is usually (slightly) faster with nstlist=0, nstype=simple and particle decomposition");
829 if (ir->nstcomm == 0)
831 ir->comm_mode = ecmNO;
833 if (ir->comm_mode != ecmNO)
837 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");
838 ir->nstcomm = abs(ir->nstcomm);
841 if (ir->nstcalcenergy > 0 && ir->nstcomm < ir->nstcalcenergy)
843 warning_note(wi, "nstcomm < nstcalcenergy defeats the purpose of nstcalcenergy, setting nstcomm to nstcalcenergy");
844 ir->nstcomm = ir->nstcalcenergy;
847 if (ir->comm_mode == ecmANGULAR)
849 sprintf(err_buf, "Can not remove the rotation around the center of mass with periodic molecules");
850 CHECK(ir->bPeriodicMols);
851 if (ir->ePBC != epbcNONE)
853 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).");
858 if (EI_STATE_VELOCITY(ir->eI) && ir->ePBC == epbcNONE && ir->comm_mode != ecmANGULAR)
860 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.");
863 sprintf(err_buf, "Twin-range neighbour searching (NS) with simple NS"
864 " algorithm not implemented");
865 CHECK(((ir->rcoulomb > ir->rlist) || (ir->rvdw > ir->rlist))
866 && (ir->ns_type == ensSIMPLE));
868 /* TEMPERATURE COUPLING */
869 if (ir->etc == etcYES)
871 ir->etc = etcBERENDSEN;
872 warning_note(wi, "Old option for temperature coupling given: "
873 "changing \"yes\" to \"Berendsen\"\n");
876 if ((ir->etc == etcNOSEHOOVER) || (ir->epc == epcMTTK))
878 if (ir->opts.nhchainlength < 1)
880 sprintf(warn_buf, "number of Nose-Hoover chains (currently %d) cannot be less than 1,reset to 1\n", ir->opts.nhchainlength);
881 ir->opts.nhchainlength = 1;
882 warning(wi, warn_buf);
885 if (ir->etc == etcNOSEHOOVER && !EI_VV(ir->eI) && ir->opts.nhchainlength > 1)
887 warning_note(wi, "leapfrog does not yet support Nose-Hoover chains, nhchainlength reset to 1");
888 ir->opts.nhchainlength = 1;
893 ir->opts.nhchainlength = 0;
896 if (ir->eI == eiVVAK)
898 sprintf(err_buf, "%s implemented primarily for validation, and requires nsttcouple = 1 and nstpcouple = 1.",
900 CHECK((ir->nsttcouple != 1) || (ir->nstpcouple != 1));
903 if (ETC_ANDERSEN(ir->etc))
905 sprintf(err_buf, "%s temperature control not supported for integrator %s.", etcoupl_names[ir->etc], ei_names[ir->eI]);
906 CHECK(!(EI_VV(ir->eI)));
908 for (i = 0; i < ir->opts.ngtc; i++)
910 sprintf(err_buf, "all tau_t must currently be equal using Andersen temperature control, violated for group %d", i);
911 CHECK(ir->opts.tau_t[0] != ir->opts.tau_t[i]);
912 sprintf(err_buf, "all tau_t must be postive using Andersen temperature control, tau_t[%d]=%10.6f",
913 i, ir->opts.tau_t[i]);
914 CHECK(ir->opts.tau_t[i] < 0);
916 if (ir->nstcomm > 0 && (ir->etc == etcANDERSEN))
918 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]);
919 warning_note(wi, warn_buf);
922 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]);
923 CHECK(ir->nstcomm > 1 && (ir->etc == etcANDERSEN));
925 for (i = 0; i < ir->opts.ngtc; i++)
927 int nsteps = (int)(ir->opts.tau_t[i]/ir->delta_t);
928 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);
929 CHECK((nsteps % ir->nstcomm) && (ir->etc == etcANDERSENMASSIVE));
932 if (ir->etc == etcBERENDSEN)
934 sprintf(warn_buf, "The %s thermostat does not generate the correct kinetic energy distribution. You might want to consider using the %s thermostat.",
935 ETCOUPLTYPE(ir->etc), ETCOUPLTYPE(etcVRESCALE));
936 warning_note(wi, warn_buf);
939 if ((ir->etc == etcNOSEHOOVER || ETC_ANDERSEN(ir->etc))
940 && ir->epc == epcBERENDSEN)
942 sprintf(warn_buf, "Using Berendsen pressure coupling invalidates the "
943 "true ensemble for the thermostat");
944 warning(wi, warn_buf);
947 /* PRESSURE COUPLING */
948 if (ir->epc == epcISOTROPIC)
950 ir->epc = epcBERENDSEN;
951 warning_note(wi, "Old option for pressure coupling given: "
952 "changing \"Isotropic\" to \"Berendsen\"\n");
955 if (ir->epc != epcNO)
957 dt_pcoupl = ir->nstpcouple*ir->delta_t;
959 sprintf(err_buf, "tau-p must be > 0 instead of %g\n", ir->tau_p);
960 CHECK(ir->tau_p <= 0);
962 if (ir->tau_p/dt_pcoupl < pcouple_min_integration_steps(ir->epc))
964 sprintf(warn_buf, "For proper integration of the %s barostat, tau-p (%g) should be at least %d times larger than nstpcouple*dt (%g)",
965 EPCOUPLTYPE(ir->epc), ir->tau_p, pcouple_min_integration_steps(ir->epc), dt_pcoupl);
966 warning(wi, warn_buf);
969 sprintf(err_buf, "compressibility must be > 0 when using pressure"
970 " coupling %s\n", EPCOUPLTYPE(ir->epc));
971 CHECK(ir->compress[XX][XX] < 0 || ir->compress[YY][YY] < 0 ||
972 ir->compress[ZZ][ZZ] < 0 ||
973 (trace(ir->compress) == 0 && ir->compress[YY][XX] <= 0 &&
974 ir->compress[ZZ][XX] <= 0 && ir->compress[ZZ][YY] <= 0));
976 if (epcPARRINELLORAHMAN == ir->epc && opts->bGenVel)
979 "You are generating velocities so I am assuming you "
980 "are equilibrating a system. You are using "
981 "%s pressure coupling, but this can be "
982 "unstable for equilibration. If your system crashes, try "
983 "equilibrating first with Berendsen pressure coupling. If "
984 "you are not equilibrating the system, you can probably "
985 "ignore this warning.",
986 epcoupl_names[ir->epc]);
987 warning(wi, warn_buf);
995 if ((ir->epc != epcBERENDSEN) && (ir->epc != epcMTTK))
997 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.");
1002 /* ELECTROSTATICS */
1003 /* More checks are in triple check (grompp.c) */
1005 if (ir->coulombtype == eelSWITCH)
1007 sprintf(warn_buf, "coulombtype = %s is only for testing purposes and can lead to serious "
1008 "artifacts, advice: use coulombtype = %s",
1009 eel_names[ir->coulombtype],
1010 eel_names[eelRF_ZERO]);
1011 warning(wi, warn_buf);
1014 if (ir->epsilon_r != 1 && ir->implicit_solvent == eisGBSA)
1016 sprintf(warn_buf, "epsilon-r = %g with GB implicit solvent, will use this value for inner dielectric", ir->epsilon_r);
1017 warning_note(wi, warn_buf);
1020 if (EEL_RF(ir->coulombtype) && ir->epsilon_rf == 1 && ir->epsilon_r != 1)
1022 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);
1023 warning(wi, warn_buf);
1024 ir->epsilon_rf = ir->epsilon_r;
1025 ir->epsilon_r = 1.0;
1028 if (getenv("GALACTIC_DYNAMICS") == NULL)
1030 sprintf(err_buf, "epsilon-r must be >= 0 instead of %g\n", ir->epsilon_r);
1031 CHECK(ir->epsilon_r < 0);
1034 if (EEL_RF(ir->coulombtype))
1036 /* reaction field (at the cut-off) */
1038 if (ir->coulombtype == eelRF_ZERO)
1040 sprintf(warn_buf, "With coulombtype = %s, epsilon-rf must be 0, assuming you meant epsilon_rf=0",
1041 eel_names[ir->coulombtype]);
1042 CHECK(ir->epsilon_rf != 0);
1043 ir->epsilon_rf = 0.0;
1046 sprintf(err_buf, "epsilon-rf must be >= epsilon-r");
1047 CHECK((ir->epsilon_rf < ir->epsilon_r && ir->epsilon_rf != 0) ||
1048 (ir->epsilon_r == 0));
1049 if (ir->epsilon_rf == ir->epsilon_r)
1051 sprintf(warn_buf, "Using epsilon-rf = epsilon-r with %s does not make sense",
1052 eel_names[ir->coulombtype]);
1053 warning(wi, warn_buf);
1056 /* Allow rlist>rcoulomb for tabulated long range stuff. This just
1057 * means the interaction is zero outside rcoulomb, but it helps to
1058 * provide accurate energy conservation.
1060 if (EEL_MIGHT_BE_ZERO_AT_CUTOFF(ir->coulombtype))
1062 if (EEL_SWITCHED(ir->coulombtype))
1065 "With coulombtype = %s rcoulomb_switch must be < rcoulomb. Or, better: Use the potential modifier options!",
1066 eel_names[ir->coulombtype]);
1067 CHECK(ir->rcoulomb_switch >= ir->rcoulomb);
1070 else if (ir->coulombtype == eelCUT || EEL_RF(ir->coulombtype))
1072 if (ir->cutoff_scheme == ecutsGROUP && ir->coulomb_modifier == eintmodNONE)
1074 sprintf(err_buf, "With coulombtype = %s, rcoulomb should be >= rlist unless you use a potential modifier",
1075 eel_names[ir->coulombtype]);
1076 CHECK(ir->rlist > ir->rcoulomb);
1080 if (ir->coulombtype == eelSWITCH || ir->coulombtype == eelSHIFT ||
1081 ir->vdwtype == evdwSWITCH || ir->vdwtype == evdwSHIFT)
1084 "The switch/shift interaction settings are just for compatibility; you will get better "
1085 "performance from applying potential modifiers to your interactions!\n");
1086 warning_note(wi, warn_buf);
1089 if (ir->coulombtype == eelPMESWITCH)
1091 if (ir->rcoulomb_switch/ir->rcoulomb < 0.9499)
1093 sprintf(warn_buf, "The switching range for %s should be 5%% or less, energy conservation will be good anyhow, since ewald_rtol = %g",
1094 eel_names[ir->coulombtype],
1096 warning(wi, warn_buf);
1100 if (EEL_FULL(ir->coulombtype))
1102 if (ir->coulombtype == eelPMESWITCH || ir->coulombtype == eelPMEUSER ||
1103 ir->coulombtype == eelPMEUSERSWITCH)
1105 sprintf(err_buf, "With coulombtype = %s, rcoulomb must be <= rlist",
1106 eel_names[ir->coulombtype]);
1107 CHECK(ir->rcoulomb > ir->rlist);
1109 else if (ir->cutoff_scheme == ecutsGROUP && ir->coulomb_modifier == eintmodNONE)
1111 if (ir->coulombtype == eelPME || ir->coulombtype == eelP3M_AD)
1114 "With coulombtype = %s (without modifier), rcoulomb must be equal to rlist,\n"
1115 "or rlistlong if nstcalclr=1. For optimal energy conservation,consider using\n"
1116 "a potential modifier.", eel_names[ir->coulombtype]);
1117 if (ir->nstcalclr == 1)
1119 CHECK(ir->rcoulomb != ir->rlist && ir->rcoulomb != ir->rlistlong);
1123 CHECK(ir->rcoulomb != ir->rlist);
1129 if (EEL_PME(ir->coulombtype))
1131 if (ir->pme_order < 3)
1133 warning_error(wi, "pme-order can not be smaller than 3");
1137 if (ir->nwall == 2 && EEL_FULL(ir->coulombtype))
1139 if (ir->ewald_geometry == eewg3D)
1141 sprintf(warn_buf, "With pbc=%s you should use ewald-geometry=%s",
1142 epbc_names[ir->ePBC], eewg_names[eewg3DC]);
1143 warning(wi, warn_buf);
1145 /* This check avoids extra pbc coding for exclusion corrections */
1146 sprintf(err_buf, "wall-ewald-zfac should be >= 2");
1147 CHECK(ir->wall_ewald_zfac < 2);
1150 if (EVDW_SWITCHED(ir->vdwtype))
1152 sprintf(err_buf, "With vdwtype = %s rvdw-switch must be < rvdw. Or, better - use a potential modifier.",
1153 evdw_names[ir->vdwtype]);
1154 CHECK(ir->rvdw_switch >= ir->rvdw);
1156 else if (ir->vdwtype == evdwCUT)
1158 if (ir->cutoff_scheme == ecutsGROUP && ir->vdw_modifier == eintmodNONE)
1160 sprintf(err_buf, "With vdwtype = %s, rvdw must be >= rlist unless you use a potential modifier", evdw_names[ir->vdwtype]);
1161 CHECK(ir->rlist > ir->rvdw);
1164 if (ir->cutoff_scheme == ecutsGROUP)
1166 if (((ir->coulomb_modifier != eintmodNONE && ir->rcoulomb == ir->rlist) ||
1167 (ir->vdw_modifier != eintmodNONE && ir->rvdw == ir->rlist)) &&
1170 warning_note(wi, "With exact cut-offs, rlist should be "
1171 "larger than rcoulomb and rvdw, so that there "
1172 "is a buffer region for particle motion "
1173 "between neighborsearch steps");
1176 if (EEL_IS_ZERO_AT_CUTOFF(ir->coulombtype)
1177 && (ir->rlistlong <= ir->rcoulomb))
1179 sprintf(warn_buf, "For energy conservation with switch/shift potentials, %s should be 0.1 to 0.3 nm larger than rcoulomb.",
1180 IR_TWINRANGE(*ir) ? "rlistlong" : "rlist");
1181 warning_note(wi, warn_buf);
1183 if (EVDW_SWITCHED(ir->vdwtype) && (ir->rlistlong <= ir->rvdw))
1185 sprintf(warn_buf, "For energy conservation with switch/shift potentials, %s should be 0.1 to 0.3 nm larger than rvdw.",
1186 IR_TWINRANGE(*ir) ? "rlistlong" : "rlist");
1187 warning_note(wi, warn_buf);
1191 if (ir->vdwtype == evdwUSER && ir->eDispCorr != edispcNO)
1193 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.");
1196 if (ir->nstlist == -1)
1198 sprintf(err_buf, "With nstlist=-1 rvdw and rcoulomb should be smaller than rlist to account for diffusion and possibly charge-group radii");
1199 CHECK(ir->rvdw >= ir->rlist || ir->rcoulomb >= ir->rlist);
1201 sprintf(err_buf, "nstlist can not be smaller than -1");
1202 CHECK(ir->nstlist < -1);
1204 if (ir->eI == eiLBFGS && (ir->coulombtype == eelCUT || ir->vdwtype == evdwCUT)
1207 warning(wi, "For efficient BFGS minimization, use switch/shift/pme instead of cut-off.");
1210 if (ir->eI == eiLBFGS && ir->nbfgscorr <= 0)
1212 warning(wi, "Using L-BFGS with nbfgscorr<=0 just gets you steepest descent.");
1215 /* ENERGY CONSERVATION */
1216 if (ir_NVE(ir) && ir->cutoff_scheme == ecutsGROUP)
1218 if (!EVDW_MIGHT_BE_ZERO_AT_CUTOFF(ir->vdwtype) && ir->rvdw > 0 && ir->vdw_modifier == eintmodNONE)
1220 sprintf(warn_buf, "You are using a cut-off for VdW interactions with NVE, for good energy conservation use vdwtype = %s (possibly with DispCorr)",
1221 evdw_names[evdwSHIFT]);
1222 warning_note(wi, warn_buf);
1224 if (!EEL_MIGHT_BE_ZERO_AT_CUTOFF(ir->coulombtype) && ir->rcoulomb > 0 && ir->coulomb_modifier == eintmodNONE)
1226 sprintf(warn_buf, "You are using a cut-off for electrostatics with NVE, for good energy conservation use coulombtype = %s or %s",
1227 eel_names[eelPMESWITCH], eel_names[eelRF_ZERO]);
1228 warning_note(wi, warn_buf);
1232 /* IMPLICIT SOLVENT */
1233 if (ir->coulombtype == eelGB_NOTUSED)
1235 ir->coulombtype = eelCUT;
1236 ir->implicit_solvent = eisGBSA;
1237 fprintf(stderr, "Note: Old option for generalized born electrostatics given:\n"
1238 "Changing coulombtype from \"generalized-born\" to \"cut-off\" and instead\n"
1239 "setting implicit-solvent value to \"GBSA\" in input section.\n");
1242 if (ir->sa_algorithm == esaSTILL)
1244 sprintf(err_buf, "Still SA algorithm not available yet, use %s or %s instead\n", esa_names[esaAPPROX], esa_names[esaNO]);
1245 CHECK(ir->sa_algorithm == esaSTILL);
1248 if (ir->implicit_solvent == eisGBSA)
1250 sprintf(err_buf, "With GBSA implicit solvent, rgbradii must be equal to rlist.");
1251 CHECK(ir->rgbradii != ir->rlist);
1253 if (ir->coulombtype != eelCUT)
1255 sprintf(err_buf, "With GBSA, coulombtype must be equal to %s\n", eel_names[eelCUT]);
1256 CHECK(ir->coulombtype != eelCUT);
1258 if (ir->vdwtype != evdwCUT)
1260 sprintf(err_buf, "With GBSA, vdw-type must be equal to %s\n", evdw_names[evdwCUT]);
1261 CHECK(ir->vdwtype != evdwCUT);
1263 if (ir->nstgbradii < 1)
1265 sprintf(warn_buf, "Using GBSA with nstgbradii<1, setting nstgbradii=1");
1266 warning_note(wi, warn_buf);
1269 if (ir->sa_algorithm == esaNO)
1271 sprintf(warn_buf, "No SA (non-polar) calculation requested together with GB. Are you sure this is what you want?\n");
1272 warning_note(wi, warn_buf);
1274 if (ir->sa_surface_tension < 0 && ir->sa_algorithm != esaNO)
1276 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");
1277 warning_note(wi, warn_buf);
1279 if (ir->gb_algorithm == egbSTILL)
1281 ir->sa_surface_tension = 0.0049 * CAL2JOULE * 100;
1285 ir->sa_surface_tension = 0.0054 * CAL2JOULE * 100;
1288 if (ir->sa_surface_tension == 0 && ir->sa_algorithm != esaNO)
1290 sprintf(err_buf, "Surface tension set to 0 while SA-calculation requested\n");
1291 CHECK(ir->sa_surface_tension == 0 && ir->sa_algorithm != esaNO);
1298 if (ir->cutoff_scheme != ecutsGROUP)
1300 warning_error(wi, "AdresS simulation supports only cutoff-scheme=group");
1304 warning_error(wi, "AdresS simulation supports only stochastic dynamics");
1306 if (ir->epc != epcNO)
1308 warning_error(wi, "AdresS simulation does not support pressure coupling");
1310 if (EEL_FULL(ir->coulombtype))
1312 warning_error(wi, "AdresS simulation does not support long-range electrostatics");
1317 /* count the number of text elemets separated by whitespace in a string.
1318 str = the input string
1319 maxptr = the maximum number of allowed elements
1320 ptr = the output array of pointers to the first character of each element
1321 returns: the number of elements. */
1322 int str_nelem(const char *str, int maxptr, char *ptr[])
1327 copy0 = strdup(str);
1330 while (*copy != '\0')
1334 gmx_fatal(FARGS, "Too many groups on line: '%s' (max is %d)",
1342 while ((*copy != '\0') && !isspace(*copy))
1361 /* interpret a number of doubles from a string and put them in an array,
1362 after allocating space for them.
1363 str = the input string
1364 n = the (pre-allocated) number of doubles read
1365 r = the output array of doubles. */
1366 static void parse_n_real(char *str, int *n, real **r)
1371 *n = str_nelem(str, MAXPTR, ptr);
1374 for (i = 0; i < *n; i++)
1376 (*r)[i] = strtod(ptr[i], NULL);
1380 static void do_fep_params(t_inputrec *ir, char fep_lambda[][STRLEN], char weights[STRLEN])
1383 int i, j, max_n_lambda, nweights, nfep[efptNR];
1384 t_lambda *fep = ir->fepvals;
1385 t_expanded *expand = ir->expandedvals;
1386 real **count_fep_lambdas;
1387 gmx_bool bOneLambda = TRUE;
1389 snew(count_fep_lambdas, efptNR);
1391 /* FEP input processing */
1392 /* first, identify the number of lambda values for each type.
1393 All that are nonzero must have the same number */
1395 for (i = 0; i < efptNR; i++)
1397 parse_n_real(fep_lambda[i], &(nfep[i]), &(count_fep_lambdas[i]));
1400 /* now, determine the number of components. All must be either zero, or equal. */
1403 for (i = 0; i < efptNR; i++)
1405 if (nfep[i] > max_n_lambda)
1407 max_n_lambda = nfep[i]; /* here's a nonzero one. All of them
1408 must have the same number if its not zero.*/
1413 for (i = 0; i < efptNR; i++)
1417 ir->fepvals->separate_dvdl[i] = FALSE;
1419 else if (nfep[i] == max_n_lambda)
1421 if (i != efptTEMPERATURE) /* we treat this differently -- not really a reason to compute the derivative with
1422 respect to the temperature currently */
1424 ir->fepvals->separate_dvdl[i] = TRUE;
1429 gmx_fatal(FARGS, "Number of lambdas (%d) for FEP type %s not equal to number of other types (%d)",
1430 nfep[i], efpt_names[i], max_n_lambda);
1433 /* we don't print out dhdl if the temperature is changing, since we can't correctly define dhdl in this case */
1434 ir->fepvals->separate_dvdl[efptTEMPERATURE] = FALSE;
1436 /* the number of lambdas is the number we've read in, which is either zero
1437 or the same for all */
1438 fep->n_lambda = max_n_lambda;
1440 /* allocate space for the array of lambda values */
1441 snew(fep->all_lambda, efptNR);
1442 /* if init_lambda is defined, we need to set lambda */
1443 if ((fep->init_lambda > 0) && (fep->n_lambda == 0))
1445 ir->fepvals->separate_dvdl[efptFEP] = TRUE;
1447 /* otherwise allocate the space for all of the lambdas, and transfer the data */
1448 for (i = 0; i < efptNR; i++)
1450 snew(fep->all_lambda[i], fep->n_lambda);
1451 if (nfep[i] > 0) /* if it's zero, then the count_fep_lambda arrays
1454 for (j = 0; j < fep->n_lambda; j++)
1456 fep->all_lambda[i][j] = (double)count_fep_lambdas[i][j];
1458 sfree(count_fep_lambdas[i]);
1461 sfree(count_fep_lambdas);
1463 /* "fep-vals" is either zero or the full number. If zero, we'll need to define fep-lambdas for internal
1464 bookkeeping -- for now, init_lambda */
1466 if ((nfep[efptFEP] == 0) && (fep->init_lambda >= 0))
1468 for (i = 0; i < fep->n_lambda; i++)
1470 fep->all_lambda[efptFEP][i] = fep->init_lambda;
1474 /* check to see if only a single component lambda is defined, and soft core is defined.
1475 In this case, turn on coulomb soft core */
1477 if (max_n_lambda == 0)
1483 for (i = 0; i < efptNR; i++)
1485 if ((nfep[i] != 0) && (i != efptFEP))
1491 if ((bOneLambda) && (fep->sc_alpha > 0))
1493 fep->bScCoul = TRUE;
1496 /* Fill in the others with the efptFEP if they are not explicitly
1497 specified (i.e. nfep[i] == 0). This means if fep is not defined,
1498 they are all zero. */
1500 for (i = 0; i < efptNR; i++)
1502 if ((nfep[i] == 0) && (i != efptFEP))
1504 for (j = 0; j < fep->n_lambda; j++)
1506 fep->all_lambda[i][j] = fep->all_lambda[efptFEP][j];
1512 /* make it easier if sc_r_power = 48 by increasing it to the 4th power, to be in the right scale. */
1513 if (fep->sc_r_power == 48)
1515 if (fep->sc_alpha > 0.1)
1517 gmx_fatal(FARGS, "sc_alpha (%f) for sc_r_power = 48 should usually be between 0.001 and 0.004", fep->sc_alpha);
1521 expand = ir->expandedvals;
1522 /* now read in the weights */
1523 parse_n_real(weights, &nweights, &(expand->init_lambda_weights));
1526 snew(expand->init_lambda_weights, fep->n_lambda); /* initialize to zero */
1528 else if (nweights != fep->n_lambda)
1530 gmx_fatal(FARGS, "Number of weights (%d) is not equal to number of lambda values (%d)",
1531 nweights, fep->n_lambda);
1533 if ((expand->nstexpanded < 0) && (ir->efep != efepNO))
1535 expand->nstexpanded = fep->nstdhdl;
1536 /* if you don't specify nstexpanded when doing expanded ensemble free energy calcs, it is set to nstdhdl */
1538 if ((expand->nstexpanded < 0) && ir->bSimTemp)
1540 expand->nstexpanded = 2*(int)(ir->opts.tau_t[0]/ir->delta_t);
1541 /* if you don't specify nstexpanded when doing expanded ensemble simulated tempering, it is set to
1542 2*tau_t just to be careful so it's not to frequent */
1547 static void do_simtemp_params(t_inputrec *ir)
1550 snew(ir->simtempvals->temperatures, ir->fepvals->n_lambda);
1551 GetSimTemps(ir->fepvals->n_lambda, ir->simtempvals, ir->fepvals->all_lambda[efptTEMPERATURE]);
1556 static void do_wall_params(t_inputrec *ir,
1557 char *wall_atomtype, char *wall_density,
1561 char *names[MAXPTR];
1564 opts->wall_atomtype[0] = NULL;
1565 opts->wall_atomtype[1] = NULL;
1567 ir->wall_atomtype[0] = -1;
1568 ir->wall_atomtype[1] = -1;
1569 ir->wall_density[0] = 0;
1570 ir->wall_density[1] = 0;
1574 nstr = str_nelem(wall_atomtype, MAXPTR, names);
1575 if (nstr != ir->nwall)
1577 gmx_fatal(FARGS, "Expected %d elements for wall_atomtype, found %d",
1580 for (i = 0; i < ir->nwall; i++)
1582 opts->wall_atomtype[i] = strdup(names[i]);
1585 if (ir->wall_type == ewt93 || ir->wall_type == ewt104)
1587 nstr = str_nelem(wall_density, MAXPTR, names);
1588 if (nstr != ir->nwall)
1590 gmx_fatal(FARGS, "Expected %d elements for wall-density, found %d", ir->nwall, nstr);
1592 for (i = 0; i < ir->nwall; i++)
1594 sscanf(names[i], "%lf", &dbl);
1597 gmx_fatal(FARGS, "wall-density[%d] = %f\n", i, dbl);
1599 ir->wall_density[i] = dbl;
1605 static void add_wall_energrps(gmx_groups_t *groups, int nwall, t_symtab *symtab)
1613 srenew(groups->grpname, groups->ngrpname+nwall);
1614 grps = &(groups->grps[egcENER]);
1615 srenew(grps->nm_ind, grps->nr+nwall);
1616 for (i = 0; i < nwall; i++)
1618 sprintf(str, "wall%d", i);
1619 groups->grpname[groups->ngrpname] = put_symtab(symtab, str);
1620 grps->nm_ind[grps->nr++] = groups->ngrpname++;
1625 void read_expandedparams(int *ninp_p, t_inpfile **inp_p,
1626 t_expanded *expand, warninp_t wi)
1628 int ninp, nerror = 0;
1634 /* read expanded ensemble parameters */
1635 CCTYPE ("expanded ensemble variables");
1636 ITYPE ("nstexpanded", expand->nstexpanded, -1);
1637 EETYPE("lmc-stats", expand->elamstats, elamstats_names);
1638 EETYPE("lmc-move", expand->elmcmove, elmcmove_names);
1639 EETYPE("lmc-weights-equil", expand->elmceq, elmceq_names);
1640 ITYPE ("weight-equil-number-all-lambda", expand->equil_n_at_lam, -1);
1641 ITYPE ("weight-equil-number-samples", expand->equil_samples, -1);
1642 ITYPE ("weight-equil-number-steps", expand->equil_steps, -1);
1643 RTYPE ("weight-equil-wl-delta", expand->equil_wl_delta, -1);
1644 RTYPE ("weight-equil-count-ratio", expand->equil_ratio, -1);
1645 CCTYPE("Seed for Monte Carlo in lambda space");
1646 ITYPE ("lmc-seed", expand->lmc_seed, -1);
1647 RTYPE ("mc-temperature", expand->mc_temp, -1);
1648 ITYPE ("lmc-repeats", expand->lmc_repeats, 1);
1649 ITYPE ("lmc-gibbsdelta", expand->gibbsdeltalam, -1);
1650 ITYPE ("lmc-forced-nstart", expand->lmc_forced_nstart, 0);
1651 EETYPE("symmetrized-transition-matrix", expand->bSymmetrizedTMatrix, yesno_names);
1652 ITYPE("nst-transition-matrix", expand->nstTij, -1);
1653 ITYPE ("mininum-var-min", expand->minvarmin, 100); /*default is reasonable */
1654 ITYPE ("weight-c-range", expand->c_range, 0); /* default is just C=0 */
1655 RTYPE ("wl-scale", expand->wl_scale, 0.8);
1656 RTYPE ("wl-ratio", expand->wl_ratio, 0.8);
1657 RTYPE ("init-wl-delta", expand->init_wl_delta, 1.0);
1658 EETYPE("wl-oneovert", expand->bWLoneovert, yesno_names);
1666 void get_ir(const char *mdparin, const char *mdparout,
1667 t_inputrec *ir, t_gromppopts *opts,
1671 double dumdub[2][6];
1675 char warn_buf[STRLEN];
1676 t_lambda *fep = ir->fepvals;
1677 t_expanded *expand = ir->expandedvals;
1679 inp = read_inpfile(mdparin, &ninp, wi);
1681 snew(dumstr[0], STRLEN);
1682 snew(dumstr[1], STRLEN);
1684 /* remove the following deprecated commands */
1687 REM_TYPE("domain-decomposition");
1688 REM_TYPE("andersen-seed");
1690 REM_TYPE("dihre-fc");
1691 REM_TYPE("dihre-tau");
1692 REM_TYPE("nstdihreout");
1693 REM_TYPE("nstcheckpoint");
1695 /* replace the following commands with the clearer new versions*/
1696 REPL_TYPE("unconstrained-start", "continuation");
1697 REPL_TYPE("foreign-lambda", "fep-lambdas");
1699 CCTYPE ("VARIOUS PREPROCESSING OPTIONS");
1700 CTYPE ("Preprocessor information: use cpp syntax.");
1701 CTYPE ("e.g.: -I/home/joe/doe -I/home/mary/roe");
1702 STYPE ("include", opts->include, NULL);
1703 CTYPE ("e.g.: -DPOSRES -DFLEXIBLE (note these variable names are case sensitive)");
1704 STYPE ("define", opts->define, NULL);
1706 CCTYPE ("RUN CONTROL PARAMETERS");
1707 EETYPE("integrator", ir->eI, ei_names);
1708 CTYPE ("Start time and timestep in ps");
1709 RTYPE ("tinit", ir->init_t, 0.0);
1710 RTYPE ("dt", ir->delta_t, 0.001);
1711 STEPTYPE ("nsteps", ir->nsteps, 0);
1712 CTYPE ("For exact run continuation or redoing part of a run");
1713 STEPTYPE ("init-step", ir->init_step, 0);
1714 CTYPE ("Part index is updated automatically on checkpointing (keeps files separate)");
1715 ITYPE ("simulation-part", ir->simulation_part, 1);
1716 CTYPE ("mode for center of mass motion removal");
1717 EETYPE("comm-mode", ir->comm_mode, ecm_names);
1718 CTYPE ("number of steps for center of mass motion removal");
1719 ITYPE ("nstcomm", ir->nstcomm, 100);
1720 CTYPE ("group(s) for center of mass motion removal");
1721 STYPE ("comm-grps", vcm, NULL);
1723 CCTYPE ("LANGEVIN DYNAMICS OPTIONS");
1724 CTYPE ("Friction coefficient (amu/ps) and random seed");
1725 RTYPE ("bd-fric", ir->bd_fric, 0.0);
1726 ITYPE ("ld-seed", ir->ld_seed, 1993);
1729 CCTYPE ("ENERGY MINIMIZATION OPTIONS");
1730 CTYPE ("Force tolerance and initial step-size");
1731 RTYPE ("emtol", ir->em_tol, 10.0);
1732 RTYPE ("emstep", ir->em_stepsize, 0.01);
1733 CTYPE ("Max number of iterations in relax-shells");
1734 ITYPE ("niter", ir->niter, 20);
1735 CTYPE ("Step size (ps^2) for minimization of flexible constraints");
1736 RTYPE ("fcstep", ir->fc_stepsize, 0);
1737 CTYPE ("Frequency of steepest descents steps when doing CG");
1738 ITYPE ("nstcgsteep", ir->nstcgsteep, 1000);
1739 ITYPE ("nbfgscorr", ir->nbfgscorr, 10);
1741 CCTYPE ("TEST PARTICLE INSERTION OPTIONS");
1742 RTYPE ("rtpi", ir->rtpi, 0.05);
1744 /* Output options */
1745 CCTYPE ("OUTPUT CONTROL OPTIONS");
1746 CTYPE ("Output frequency for coords (x), velocities (v) and forces (f)");
1747 ITYPE ("nstxout", ir->nstxout, 0);
1748 ITYPE ("nstvout", ir->nstvout, 0);
1749 ITYPE ("nstfout", ir->nstfout, 0);
1750 ir->nstcheckpoint = 1000;
1751 CTYPE ("Output frequency for energies to log file and energy file");
1752 ITYPE ("nstlog", ir->nstlog, 1000);
1753 ITYPE ("nstcalcenergy", ir->nstcalcenergy, 100);
1754 ITYPE ("nstenergy", ir->nstenergy, 1000);
1755 CTYPE ("Output frequency and precision for .xtc file");
1756 ITYPE ("nstxtcout", ir->nstxtcout, 0);
1757 RTYPE ("xtc-precision", ir->xtcprec, 1000.0);
1758 CTYPE ("This selects the subset of atoms for the .xtc file. You can");
1759 CTYPE ("select multiple groups. By default all atoms will be written.");
1760 STYPE ("xtc-grps", xtc_grps, NULL);
1761 CTYPE ("Selection of energy groups");
1762 STYPE ("energygrps", energy, NULL);
1764 /* Neighbor searching */
1765 CCTYPE ("NEIGHBORSEARCHING PARAMETERS");
1766 CTYPE ("cut-off scheme (group: using charge groups, Verlet: particle based cut-offs)");
1767 EETYPE("cutoff-scheme", ir->cutoff_scheme, ecutscheme_names);
1768 CTYPE ("nblist update frequency");
1769 ITYPE ("nstlist", ir->nstlist, 10);
1770 CTYPE ("ns algorithm (simple or grid)");
1771 EETYPE("ns-type", ir->ns_type, ens_names);
1772 /* set ndelta to the optimal value of 2 */
1774 CTYPE ("Periodic boundary conditions: xyz, no, xy");
1775 EETYPE("pbc", ir->ePBC, epbc_names);
1776 EETYPE("periodic-molecules", ir->bPeriodicMols, yesno_names);
1777 CTYPE ("Allowed energy drift due to the Verlet buffer in kJ/mol/ps per atom,");
1778 CTYPE ("a value of -1 means: use rlist");
1779 RTYPE("verlet-buffer-drift", ir->verletbuf_drift, 0.005);
1780 CTYPE ("nblist cut-off");
1781 RTYPE ("rlist", ir->rlist, 1.0);
1782 CTYPE ("long-range cut-off for switched potentials");
1783 RTYPE ("rlistlong", ir->rlistlong, -1);
1784 ITYPE ("nstcalclr", ir->nstcalclr, -1);
1786 /* Electrostatics */
1787 CCTYPE ("OPTIONS FOR ELECTROSTATICS AND VDW");
1788 CTYPE ("Method for doing electrostatics");
1789 EETYPE("coulombtype", ir->coulombtype, eel_names);
1790 EETYPE("coulomb-modifier", ir->coulomb_modifier, eintmod_names);
1791 CTYPE ("cut-off lengths");
1792 RTYPE ("rcoulomb-switch", ir->rcoulomb_switch, 0.0);
1793 RTYPE ("rcoulomb", ir->rcoulomb, 1.0);
1794 CTYPE ("Relative dielectric constant for the medium and the reaction field");
1795 RTYPE ("epsilon-r", ir->epsilon_r, 1.0);
1796 RTYPE ("epsilon-rf", ir->epsilon_rf, 0.0);
1797 CTYPE ("Method for doing Van der Waals");
1798 EETYPE("vdw-type", ir->vdwtype, evdw_names);
1799 EETYPE("vdw-modifier", ir->vdw_modifier, eintmod_names);
1800 CTYPE ("cut-off lengths");
1801 RTYPE ("rvdw-switch", ir->rvdw_switch, 0.0);
1802 RTYPE ("rvdw", ir->rvdw, 1.0);
1803 CTYPE ("Apply long range dispersion corrections for Energy and Pressure");
1804 EETYPE("DispCorr", ir->eDispCorr, edispc_names);
1805 CTYPE ("Extension of the potential lookup tables beyond the cut-off");
1806 RTYPE ("table-extension", ir->tabext, 1.0);
1807 CTYPE ("Separate tables between energy group pairs");
1808 STYPE ("energygrp-table", egptable, NULL);
1809 CTYPE ("Spacing for the PME/PPPM FFT grid");
1810 RTYPE ("fourierspacing", ir->fourier_spacing, 0.12);
1811 CTYPE ("FFT grid size, when a value is 0 fourierspacing will be used");
1812 ITYPE ("fourier-nx", ir->nkx, 0);
1813 ITYPE ("fourier-ny", ir->nky, 0);
1814 ITYPE ("fourier-nz", ir->nkz, 0);
1815 CTYPE ("EWALD/PME/PPPM parameters");
1816 ITYPE ("pme-order", ir->pme_order, 4);
1817 RTYPE ("ewald-rtol", ir->ewald_rtol, 0.00001);
1818 EETYPE("ewald-geometry", ir->ewald_geometry, eewg_names);
1819 RTYPE ("epsilon-surface", ir->epsilon_surface, 0.0);
1820 EETYPE("optimize-fft", ir->bOptFFT, yesno_names);
1822 CCTYPE("IMPLICIT SOLVENT ALGORITHM");
1823 EETYPE("implicit-solvent", ir->implicit_solvent, eis_names);
1825 CCTYPE ("GENERALIZED BORN ELECTROSTATICS");
1826 CTYPE ("Algorithm for calculating Born radii");
1827 EETYPE("gb-algorithm", ir->gb_algorithm, egb_names);
1828 CTYPE ("Frequency of calculating the Born radii inside rlist");
1829 ITYPE ("nstgbradii", ir->nstgbradii, 1);
1830 CTYPE ("Cutoff for Born radii calculation; the contribution from atoms");
1831 CTYPE ("between rlist and rgbradii is updated every nstlist steps");
1832 RTYPE ("rgbradii", ir->rgbradii, 1.0);
1833 CTYPE ("Dielectric coefficient of the implicit solvent");
1834 RTYPE ("gb-epsilon-solvent", ir->gb_epsilon_solvent, 80.0);
1835 CTYPE ("Salt concentration in M for Generalized Born models");
1836 RTYPE ("gb-saltconc", ir->gb_saltconc, 0.0);
1837 CTYPE ("Scaling factors used in the OBC GB model. Default values are OBC(II)");
1838 RTYPE ("gb-obc-alpha", ir->gb_obc_alpha, 1.0);
1839 RTYPE ("gb-obc-beta", ir->gb_obc_beta, 0.8);
1840 RTYPE ("gb-obc-gamma", ir->gb_obc_gamma, 4.85);
1841 RTYPE ("gb-dielectric-offset", ir->gb_dielectric_offset, 0.009);
1842 EETYPE("sa-algorithm", ir->sa_algorithm, esa_names);
1843 CTYPE ("Surface tension (kJ/mol/nm^2) for the SA (nonpolar surface) part of GBSA");
1844 CTYPE ("The value -1 will set default value for Still/HCT/OBC GB-models.");
1845 RTYPE ("sa-surface-tension", ir->sa_surface_tension, -1);
1847 /* Coupling stuff */
1848 CCTYPE ("OPTIONS FOR WEAK COUPLING ALGORITHMS");
1849 CTYPE ("Temperature coupling");
1850 EETYPE("tcoupl", ir->etc, etcoupl_names);
1851 ITYPE ("nsttcouple", ir->nsttcouple, -1);
1852 ITYPE("nh-chain-length", ir->opts.nhchainlength, NHCHAINLENGTH);
1853 EETYPE("print-nose-hoover-chain-variables", ir->bPrintNHChains, yesno_names);
1854 CTYPE ("Groups to couple separately");
1855 STYPE ("tc-grps", tcgrps, NULL);
1856 CTYPE ("Time constant (ps) and reference temperature (K)");
1857 STYPE ("tau-t", tau_t, NULL);
1858 STYPE ("ref-t", ref_t, NULL);
1859 CTYPE ("pressure coupling");
1860 EETYPE("pcoupl", ir->epc, epcoupl_names);
1861 EETYPE("pcoupltype", ir->epct, epcoupltype_names);
1862 ITYPE ("nstpcouple", ir->nstpcouple, -1);
1863 CTYPE ("Time constant (ps), compressibility (1/bar) and reference P (bar)");
1864 RTYPE ("tau-p", ir->tau_p, 1.0);
1865 STYPE ("compressibility", dumstr[0], NULL);
1866 STYPE ("ref-p", dumstr[1], NULL);
1867 CTYPE ("Scaling of reference coordinates, No, All or COM");
1868 EETYPE ("refcoord-scaling", ir->refcoord_scaling, erefscaling_names);
1871 CCTYPE ("OPTIONS FOR QMMM calculations");
1872 EETYPE("QMMM", ir->bQMMM, yesno_names);
1873 CTYPE ("Groups treated Quantum Mechanically");
1874 STYPE ("QMMM-grps", QMMM, NULL);
1875 CTYPE ("QM method");
1876 STYPE("QMmethod", QMmethod, NULL);
1877 CTYPE ("QMMM scheme");
1878 EETYPE("QMMMscheme", ir->QMMMscheme, eQMMMscheme_names);
1879 CTYPE ("QM basisset");
1880 STYPE("QMbasis", QMbasis, NULL);
1881 CTYPE ("QM charge");
1882 STYPE ("QMcharge", QMcharge, NULL);
1883 CTYPE ("QM multiplicity");
1884 STYPE ("QMmult", QMmult, NULL);
1885 CTYPE ("Surface Hopping");
1886 STYPE ("SH", bSH, NULL);
1887 CTYPE ("CAS space options");
1888 STYPE ("CASorbitals", CASorbitals, NULL);
1889 STYPE ("CASelectrons", CASelectrons, NULL);
1890 STYPE ("SAon", SAon, NULL);
1891 STYPE ("SAoff", SAoff, NULL);
1892 STYPE ("SAsteps", SAsteps, NULL);
1893 CTYPE ("Scale factor for MM charges");
1894 RTYPE ("MMChargeScaleFactor", ir->scalefactor, 1.0);
1895 CTYPE ("Optimization of QM subsystem");
1896 STYPE ("bOPT", bOPT, NULL);
1897 STYPE ("bTS", bTS, NULL);
1899 /* Simulated annealing */
1900 CCTYPE("SIMULATED ANNEALING");
1901 CTYPE ("Type of annealing for each temperature group (no/single/periodic)");
1902 STYPE ("annealing", anneal, NULL);
1903 CTYPE ("Number of time points to use for specifying annealing in each group");
1904 STYPE ("annealing-npoints", anneal_npoints, NULL);
1905 CTYPE ("List of times at the annealing points for each group");
1906 STYPE ("annealing-time", anneal_time, NULL);
1907 CTYPE ("Temp. at each annealing point, for each group.");
1908 STYPE ("annealing-temp", anneal_temp, NULL);
1911 CCTYPE ("GENERATE VELOCITIES FOR STARTUP RUN");
1912 EETYPE("gen-vel", opts->bGenVel, yesno_names);
1913 RTYPE ("gen-temp", opts->tempi, 300.0);
1914 ITYPE ("gen-seed", opts->seed, 173529);
1917 CCTYPE ("OPTIONS FOR BONDS");
1918 EETYPE("constraints", opts->nshake, constraints);
1919 CTYPE ("Type of constraint algorithm");
1920 EETYPE("constraint-algorithm", ir->eConstrAlg, econstr_names);
1921 CTYPE ("Do not constrain the start configuration");
1922 EETYPE("continuation", ir->bContinuation, yesno_names);
1923 CTYPE ("Use successive overrelaxation to reduce the number of shake iterations");
1924 EETYPE("Shake-SOR", ir->bShakeSOR, yesno_names);
1925 CTYPE ("Relative tolerance of shake");
1926 RTYPE ("shake-tol", ir->shake_tol, 0.0001);
1927 CTYPE ("Highest order in the expansion of the constraint coupling matrix");
1928 ITYPE ("lincs-order", ir->nProjOrder, 4);
1929 CTYPE ("Number of iterations in the final step of LINCS. 1 is fine for");
1930 CTYPE ("normal simulations, but use 2 to conserve energy in NVE runs.");
1931 CTYPE ("For energy minimization with constraints it should be 4 to 8.");
1932 ITYPE ("lincs-iter", ir->nLincsIter, 1);
1933 CTYPE ("Lincs will write a warning to the stderr if in one step a bond");
1934 CTYPE ("rotates over more degrees than");
1935 RTYPE ("lincs-warnangle", ir->LincsWarnAngle, 30.0);
1936 CTYPE ("Convert harmonic bonds to morse potentials");
1937 EETYPE("morse", opts->bMorse, yesno_names);
1939 /* Energy group exclusions */
1940 CCTYPE ("ENERGY GROUP EXCLUSIONS");
1941 CTYPE ("Pairs of energy groups for which all non-bonded interactions are excluded");
1942 STYPE ("energygrp-excl", egpexcl, NULL);
1946 CTYPE ("Number of walls, type, atom types, densities and box-z scale factor for Ewald");
1947 ITYPE ("nwall", ir->nwall, 0);
1948 EETYPE("wall-type", ir->wall_type, ewt_names);
1949 RTYPE ("wall-r-linpot", ir->wall_r_linpot, -1);
1950 STYPE ("wall-atomtype", wall_atomtype, NULL);
1951 STYPE ("wall-density", wall_density, NULL);
1952 RTYPE ("wall-ewald-zfac", ir->wall_ewald_zfac, 3);
1955 CCTYPE("COM PULLING");
1956 CTYPE("Pull type: no, umbrella, constraint or constant-force");
1957 EETYPE("pull", ir->ePull, epull_names);
1958 if (ir->ePull != epullNO)
1961 pull_grp = read_pullparams(&ninp, &inp, ir->pull, &opts->pull_start, wi);
1964 /* Enforced rotation */
1965 CCTYPE("ENFORCED ROTATION");
1966 CTYPE("Enforced rotation: No or Yes");
1967 EETYPE("rotation", ir->bRot, yesno_names);
1971 rot_grp = read_rotparams(&ninp, &inp, ir->rot, wi);
1975 CCTYPE("NMR refinement stuff");
1976 CTYPE ("Distance restraints type: No, Simple or Ensemble");
1977 EETYPE("disre", ir->eDisre, edisre_names);
1978 CTYPE ("Force weighting of pairs in one distance restraint: Conservative or Equal");
1979 EETYPE("disre-weighting", ir->eDisreWeighting, edisreweighting_names);
1980 CTYPE ("Use sqrt of the time averaged times the instantaneous violation");
1981 EETYPE("disre-mixed", ir->bDisreMixed, yesno_names);
1982 RTYPE ("disre-fc", ir->dr_fc, 1000.0);
1983 RTYPE ("disre-tau", ir->dr_tau, 0.0);
1984 CTYPE ("Output frequency for pair distances to energy file");
1985 ITYPE ("nstdisreout", ir->nstdisreout, 100);
1986 CTYPE ("Orientation restraints: No or Yes");
1987 EETYPE("orire", opts->bOrire, yesno_names);
1988 CTYPE ("Orientation restraints force constant and tau for time averaging");
1989 RTYPE ("orire-fc", ir->orires_fc, 0.0);
1990 RTYPE ("orire-tau", ir->orires_tau, 0.0);
1991 STYPE ("orire-fitgrp", orirefitgrp, NULL);
1992 CTYPE ("Output frequency for trace(SD) and S to energy file");
1993 ITYPE ("nstorireout", ir->nstorireout, 100);
1995 /* free energy variables */
1996 CCTYPE ("Free energy variables");
1997 EETYPE("free-energy", ir->efep, efep_names);
1998 STYPE ("couple-moltype", couple_moltype, NULL);
1999 EETYPE("couple-lambda0", opts->couple_lam0, couple_lam);
2000 EETYPE("couple-lambda1", opts->couple_lam1, couple_lam);
2001 EETYPE("couple-intramol", opts->bCoupleIntra, yesno_names);
2003 RTYPE ("init-lambda", fep->init_lambda, -1); /* start with -1 so
2005 it was not entered */
2006 ITYPE ("init-lambda-state", fep->init_fep_state, -1);
2007 RTYPE ("delta-lambda", fep->delta_lambda, 0.0);
2008 ITYPE ("nstdhdl", fep->nstdhdl, 50);
2009 STYPE ("fep-lambdas", fep_lambda[efptFEP], NULL);
2010 STYPE ("mass-lambdas", fep_lambda[efptMASS], NULL);
2011 STYPE ("coul-lambdas", fep_lambda[efptCOUL], NULL);
2012 STYPE ("vdw-lambdas", fep_lambda[efptVDW], NULL);
2013 STYPE ("bonded-lambdas", fep_lambda[efptBONDED], NULL);
2014 STYPE ("restraint-lambdas", fep_lambda[efptRESTRAINT], NULL);
2015 STYPE ("temperature-lambdas", fep_lambda[efptTEMPERATURE], NULL);
2016 ITYPE ("calc-lambda-neighbors", fep->lambda_neighbors, 1);
2017 STYPE ("init-lambda-weights", lambda_weights, NULL);
2018 EETYPE("dhdl-print-energy", fep->bPrintEnergy, yesno_names);
2019 RTYPE ("sc-alpha", fep->sc_alpha, 0.0);
2020 ITYPE ("sc-power", fep->sc_power, 1);
2021 RTYPE ("sc-r-power", fep->sc_r_power, 6.0);
2022 RTYPE ("sc-sigma", fep->sc_sigma, 0.3);
2023 EETYPE("sc-coul", fep->bScCoul, yesno_names);
2024 ITYPE ("dh_hist_size", fep->dh_hist_size, 0);
2025 RTYPE ("dh_hist_spacing", fep->dh_hist_spacing, 0.1);
2026 EETYPE("separate-dhdl-file", fep->separate_dhdl_file,
2027 separate_dhdl_file_names);
2028 EETYPE("dhdl-derivatives", fep->dhdl_derivatives, dhdl_derivatives_names);
2029 ITYPE ("dh_hist_size", fep->dh_hist_size, 0);
2030 RTYPE ("dh_hist_spacing", fep->dh_hist_spacing, 0.1);
2032 /* Non-equilibrium MD stuff */
2033 CCTYPE("Non-equilibrium MD stuff");
2034 STYPE ("acc-grps", accgrps, NULL);
2035 STYPE ("accelerate", acc, NULL);
2036 STYPE ("freezegrps", freeze, NULL);
2037 STYPE ("freezedim", frdim, NULL);
2038 RTYPE ("cos-acceleration", ir->cos_accel, 0);
2039 STYPE ("deform", deform, NULL);
2041 /* simulated tempering variables */
2042 CCTYPE("simulated tempering variables");
2043 EETYPE("simulated-tempering", ir->bSimTemp, yesno_names);
2044 EETYPE("simulated-tempering-scaling", ir->simtempvals->eSimTempScale, esimtemp_names);
2045 RTYPE("sim-temp-low", ir->simtempvals->simtemp_low, 300.0);
2046 RTYPE("sim-temp-high", ir->simtempvals->simtemp_high, 300.0);
2048 /* expanded ensemble variables */
2049 if (ir->efep == efepEXPANDED || ir->bSimTemp)
2051 read_expandedparams(&ninp, &inp, expand, wi);
2054 /* Electric fields */
2055 CCTYPE("Electric fields");
2056 CTYPE ("Format is number of terms (int) and for all terms an amplitude (real)");
2057 CTYPE ("and a phase angle (real)");
2058 STYPE ("E-x", efield_x, NULL);
2059 STYPE ("E-xt", efield_xt, NULL);
2060 STYPE ("E-y", efield_y, NULL);
2061 STYPE ("E-yt", efield_yt, NULL);
2062 STYPE ("E-z", efield_z, NULL);
2063 STYPE ("E-zt", efield_zt, NULL);
2065 /* AdResS defined thingies */
2066 CCTYPE ("AdResS parameters");
2067 EETYPE("adress", ir->bAdress, yesno_names);
2070 snew(ir->adress, 1);
2071 read_adressparams(&ninp, &inp, ir->adress, wi);
2074 /* User defined thingies */
2075 CCTYPE ("User defined thingies");
2076 STYPE ("user1-grps", user1, NULL);
2077 STYPE ("user2-grps", user2, NULL);
2078 ITYPE ("userint1", ir->userint1, 0);
2079 ITYPE ("userint2", ir->userint2, 0);
2080 ITYPE ("userint3", ir->userint3, 0);
2081 ITYPE ("userint4", ir->userint4, 0);
2082 RTYPE ("userreal1", ir->userreal1, 0);
2083 RTYPE ("userreal2", ir->userreal2, 0);
2084 RTYPE ("userreal3", ir->userreal3, 0);
2085 RTYPE ("userreal4", ir->userreal4, 0);
2088 write_inpfile(mdparout, ninp, inp, FALSE, wi);
2089 for (i = 0; (i < ninp); i++)
2092 sfree(inp[i].value);
2096 /* Process options if necessary */
2097 for (m = 0; m < 2; m++)
2099 for (i = 0; i < 2*DIM; i++)
2108 if (sscanf(dumstr[m], "%lf", &(dumdub[m][XX])) != 1)
2110 warning_error(wi, "Pressure coupling not enough values (I need 1)");
2112 dumdub[m][YY] = dumdub[m][ZZ] = dumdub[m][XX];
2114 case epctSEMIISOTROPIC:
2115 case epctSURFACETENSION:
2116 if (sscanf(dumstr[m], "%lf%lf",
2117 &(dumdub[m][XX]), &(dumdub[m][ZZ])) != 2)
2119 warning_error(wi, "Pressure coupling not enough values (I need 2)");
2121 dumdub[m][YY] = dumdub[m][XX];
2123 case epctANISOTROPIC:
2124 if (sscanf(dumstr[m], "%lf%lf%lf%lf%lf%lf",
2125 &(dumdub[m][XX]), &(dumdub[m][YY]), &(dumdub[m][ZZ]),
2126 &(dumdub[m][3]), &(dumdub[m][4]), &(dumdub[m][5])) != 6)
2128 warning_error(wi, "Pressure coupling not enough values (I need 6)");
2132 gmx_fatal(FARGS, "Pressure coupling type %s not implemented yet",
2133 epcoupltype_names[ir->epct]);
2137 clear_mat(ir->ref_p);
2138 clear_mat(ir->compress);
2139 for (i = 0; i < DIM; i++)
2141 ir->ref_p[i][i] = dumdub[1][i];
2142 ir->compress[i][i] = dumdub[0][i];
2144 if (ir->epct == epctANISOTROPIC)
2146 ir->ref_p[XX][YY] = dumdub[1][3];
2147 ir->ref_p[XX][ZZ] = dumdub[1][4];
2148 ir->ref_p[YY][ZZ] = dumdub[1][5];
2149 if (ir->ref_p[XX][YY] != 0 && ir->ref_p[XX][ZZ] != 0 && ir->ref_p[YY][ZZ] != 0)
2151 warning(wi, "All off-diagonal reference pressures are non-zero. Are you sure you want to apply a threefold shear stress?\n");
2153 ir->compress[XX][YY] = dumdub[0][3];
2154 ir->compress[XX][ZZ] = dumdub[0][4];
2155 ir->compress[YY][ZZ] = dumdub[0][5];
2156 for (i = 0; i < DIM; i++)
2158 for (m = 0; m < i; m++)
2160 ir->ref_p[i][m] = ir->ref_p[m][i];
2161 ir->compress[i][m] = ir->compress[m][i];
2166 if (ir->comm_mode == ecmNO)
2171 opts->couple_moltype = NULL;
2172 if (strlen(couple_moltype) > 0)
2174 if (ir->efep != efepNO)
2176 opts->couple_moltype = strdup(couple_moltype);
2177 if (opts->couple_lam0 == opts->couple_lam1)
2179 warning(wi, "The lambda=0 and lambda=1 states for coupling are identical");
2181 if (ir->eI == eiMD && (opts->couple_lam0 == ecouplamNONE ||
2182 opts->couple_lam1 == ecouplamNONE))
2184 warning(wi, "For proper sampling of the (nearly) decoupled state, stochastic dynamics should be used");
2189 warning(wi, "Can not couple a molecule with free_energy = no");
2192 /* FREE ENERGY AND EXPANDED ENSEMBLE OPTIONS */
2193 if (ir->efep != efepNO)
2195 if (fep->delta_lambda > 0)
2197 ir->efep = efepSLOWGROWTH;
2203 fep->bPrintEnergy = TRUE;
2204 /* always print out the energy to dhdl if we are doing expanded ensemble, since we need the total energy
2205 if the temperature is changing. */
2208 if ((ir->efep != efepNO) || ir->bSimTemp)
2210 ir->bExpanded = FALSE;
2211 if ((ir->efep == efepEXPANDED) || ir->bSimTemp)
2213 ir->bExpanded = TRUE;
2215 do_fep_params(ir, fep_lambda, lambda_weights);
2216 if (ir->bSimTemp) /* done after fep params */
2218 do_simtemp_params(ir);
2223 ir->fepvals->n_lambda = 0;
2226 /* WALL PARAMETERS */
2228 do_wall_params(ir, wall_atomtype, wall_density, opts);
2230 /* ORIENTATION RESTRAINT PARAMETERS */
2232 if (opts->bOrire && str_nelem(orirefitgrp, MAXPTR, NULL) != 1)
2234 warning_error(wi, "ERROR: Need one orientation restraint fit group\n");
2237 /* DEFORMATION PARAMETERS */
2239 clear_mat(ir->deform);
2240 for (i = 0; i < 6; i++)
2244 m = sscanf(deform, "%lf %lf %lf %lf %lf %lf",
2245 &(dumdub[0][0]), &(dumdub[0][1]), &(dumdub[0][2]),
2246 &(dumdub[0][3]), &(dumdub[0][4]), &(dumdub[0][5]));
2247 for (i = 0; i < 3; i++)
2249 ir->deform[i][i] = dumdub[0][i];
2251 ir->deform[YY][XX] = dumdub[0][3];
2252 ir->deform[ZZ][XX] = dumdub[0][4];
2253 ir->deform[ZZ][YY] = dumdub[0][5];
2254 if (ir->epc != epcNO)
2256 for (i = 0; i < 3; i++)
2258 for (j = 0; j <= i; j++)
2260 if (ir->deform[i][j] != 0 && ir->compress[i][j] != 0)
2262 warning_error(wi, "A box element has deform set and compressibility > 0");
2266 for (i = 0; i < 3; i++)
2268 for (j = 0; j < i; j++)
2270 if (ir->deform[i][j] != 0)
2272 for (m = j; m < DIM; m++)
2274 if (ir->compress[m][j] != 0)
2276 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.");
2277 warning(wi, warn_buf);
2289 static int search_QMstring(char *s, int ng, const char *gn[])
2291 /* same as normal search_string, but this one searches QM strings */
2294 for (i = 0; (i < ng); i++)
2296 if (gmx_strcasecmp(s, gn[i]) == 0)
2302 gmx_fatal(FARGS, "this QM method or basisset (%s) is not implemented\n!", s);
2306 } /* search_QMstring */
2309 int search_string(char *s, int ng, char *gn[])
2313 for (i = 0; (i < ng); i++)
2315 if (gmx_strcasecmp(s, gn[i]) == 0)
2322 "Group %s referenced in the .mdp file was not found in the index file.\n"
2323 "Group names must match either [moleculetype] names or custom index group\n"
2324 "names, in which case you must supply an index file to the '-n' option\n"
2331 static gmx_bool do_numbering(int natoms, gmx_groups_t *groups, int ng, char *ptrs[],
2332 t_blocka *block, char *gnames[],
2333 int gtype, int restnm,
2334 int grptp, gmx_bool bVerbose,
2337 unsigned short *cbuf;
2338 t_grps *grps = &(groups->grps[gtype]);
2339 int i, j, gid, aj, ognr, ntot = 0;
2342 char warn_buf[STRLEN];
2346 fprintf(debug, "Starting numbering %d groups of type %d\n", ng, gtype);
2349 title = gtypes[gtype];
2352 /* Mark all id's as not set */
2353 for (i = 0; (i < natoms); i++)
2358 snew(grps->nm_ind, ng+1); /* +1 for possible rest group */
2359 for (i = 0; (i < ng); i++)
2361 /* Lookup the group name in the block structure */
2362 gid = search_string(ptrs[i], block->nr, gnames);
2363 if ((grptp != egrptpONE) || (i == 0))
2365 grps->nm_ind[grps->nr++] = gid;
2369 fprintf(debug, "Found gid %d for group %s\n", gid, ptrs[i]);
2372 /* Now go over the atoms in the group */
2373 for (j = block->index[gid]; (j < block->index[gid+1]); j++)
2378 /* Range checking */
2379 if ((aj < 0) || (aj >= natoms))
2381 gmx_fatal(FARGS, "Invalid atom number %d in indexfile", aj);
2383 /* Lookup up the old group number */
2387 gmx_fatal(FARGS, "Atom %d in multiple %s groups (%d and %d)",
2388 aj+1, title, ognr+1, i+1);
2392 /* Store the group number in buffer */
2393 if (grptp == egrptpONE)
2406 /* Now check whether we have done all atoms */
2410 if (grptp == egrptpALL)
2412 gmx_fatal(FARGS, "%d atoms are not part of any of the %s groups",
2413 natoms-ntot, title);
2415 else if (grptp == egrptpPART)
2417 sprintf(warn_buf, "%d atoms are not part of any of the %s groups",
2418 natoms-ntot, title);
2419 warning_note(wi, warn_buf);
2421 /* Assign all atoms currently unassigned to a rest group */
2422 for (j = 0; (j < natoms); j++)
2424 if (cbuf[j] == NOGID)
2430 if (grptp != egrptpPART)
2435 "Making dummy/rest group for %s containing %d elements\n",
2436 title, natoms-ntot);
2438 /* Add group name "rest" */
2439 grps->nm_ind[grps->nr] = restnm;
2441 /* Assign the rest name to all atoms not currently assigned to a group */
2442 for (j = 0; (j < natoms); j++)
2444 if (cbuf[j] == NOGID)
2453 if (grps->nr == 1 && (ntot == 0 || ntot == natoms))
2455 /* All atoms are part of one (or no) group, no index required */
2456 groups->ngrpnr[gtype] = 0;
2457 groups->grpnr[gtype] = NULL;
2461 groups->ngrpnr[gtype] = natoms;
2462 snew(groups->grpnr[gtype], natoms);
2463 for (j = 0; (j < natoms); j++)
2465 groups->grpnr[gtype][j] = cbuf[j];
2471 return (bRest && grptp == egrptpPART);
2474 static void calc_nrdf(gmx_mtop_t *mtop, t_inputrec *ir, char **gnames)
2477 gmx_groups_t *groups;
2479 int natoms, ai, aj, i, j, d, g, imin, jmin, nc;
2481 int *nrdf2, *na_vcm, na_tot;
2482 double *nrdf_tc, *nrdf_vcm, nrdf_uc, n_sub = 0;
2483 gmx_mtop_atomloop_all_t aloop;
2485 int mb, mol, ftype, as;
2486 gmx_molblock_t *molb;
2487 gmx_moltype_t *molt;
2490 * First calc 3xnr-atoms for each group
2491 * then subtract half a degree of freedom for each constraint
2493 * Only atoms and nuclei contribute to the degrees of freedom...
2498 groups = &mtop->groups;
2499 natoms = mtop->natoms;
2501 /* Allocate one more for a possible rest group */
2502 /* We need to sum degrees of freedom into doubles,
2503 * since floats give too low nrdf's above 3 million atoms.
2505 snew(nrdf_tc, groups->grps[egcTC].nr+1);
2506 snew(nrdf_vcm, groups->grps[egcVCM].nr+1);
2507 snew(na_vcm, groups->grps[egcVCM].nr+1);
2509 for (i = 0; i < groups->grps[egcTC].nr; i++)
2513 for (i = 0; i < groups->grps[egcVCM].nr+1; i++)
2518 snew(nrdf2, natoms);
2519 aloop = gmx_mtop_atomloop_all_init(mtop);
2520 while (gmx_mtop_atomloop_all_next(aloop, &i, &atom))
2523 if (atom->ptype == eptAtom || atom->ptype == eptNucleus)
2525 g = ggrpnr(groups, egcFREEZE, i);
2526 /* Double count nrdf for particle i */
2527 for (d = 0; d < DIM; d++)
2529 if (opts->nFreeze[g][d] == 0)
2534 nrdf_tc [ggrpnr(groups, egcTC, i)] += 0.5*nrdf2[i];
2535 nrdf_vcm[ggrpnr(groups, egcVCM, i)] += 0.5*nrdf2[i];
2540 for (mb = 0; mb < mtop->nmolblock; mb++)
2542 molb = &mtop->molblock[mb];
2543 molt = &mtop->moltype[molb->type];
2544 atom = molt->atoms.atom;
2545 for (mol = 0; mol < molb->nmol; mol++)
2547 for (ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++)
2549 ia = molt->ilist[ftype].iatoms;
2550 for (i = 0; i < molt->ilist[ftype].nr; )
2552 /* Subtract degrees of freedom for the constraints,
2553 * if the particles still have degrees of freedom left.
2554 * If one of the particles is a vsite or a shell, then all
2555 * constraint motion will go there, but since they do not
2556 * contribute to the constraints the degrees of freedom do not
2561 if (((atom[ia[1]].ptype == eptNucleus) ||
2562 (atom[ia[1]].ptype == eptAtom)) &&
2563 ((atom[ia[2]].ptype == eptNucleus) ||
2564 (atom[ia[2]].ptype == eptAtom)))
2582 imin = min(imin, nrdf2[ai]);
2583 jmin = min(jmin, nrdf2[aj]);
2586 nrdf_tc [ggrpnr(groups, egcTC, ai)] -= 0.5*imin;
2587 nrdf_tc [ggrpnr(groups, egcTC, aj)] -= 0.5*jmin;
2588 nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5*imin;
2589 nrdf_vcm[ggrpnr(groups, egcVCM, aj)] -= 0.5*jmin;
2591 ia += interaction_function[ftype].nratoms+1;
2592 i += interaction_function[ftype].nratoms+1;
2595 ia = molt->ilist[F_SETTLE].iatoms;
2596 for (i = 0; i < molt->ilist[F_SETTLE].nr; )
2598 /* Subtract 1 dof from every atom in the SETTLE */
2599 for (j = 0; j < 3; j++)
2602 imin = min(2, nrdf2[ai]);
2604 nrdf_tc [ggrpnr(groups, egcTC, ai)] -= 0.5*imin;
2605 nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5*imin;
2610 as += molt->atoms.nr;
2614 if (ir->ePull == epullCONSTRAINT)
2616 /* Correct nrdf for the COM constraints.
2617 * We correct using the TC and VCM group of the first atom
2618 * in the reference and pull group. If atoms in one pull group
2619 * belong to different TC or VCM groups it is anyhow difficult
2620 * to determine the optimal nrdf assignment.
2623 if (pull->eGeom == epullgPOS)
2626 for (i = 0; i < DIM; i++)
2638 for (i = 0; i < pull->ngrp; i++)
2641 if (pull->grp[0].nat > 0)
2643 /* Subtract 1/2 dof from the reference group */
2644 ai = pull->grp[0].ind[0];
2645 if (nrdf_tc[ggrpnr(groups, egcTC, ai)] > 1)
2647 nrdf_tc [ggrpnr(groups, egcTC, ai)] -= 0.5;
2648 nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5;
2652 /* Subtract 1/2 dof from the pulled group */
2653 ai = pull->grp[1+i].ind[0];
2654 nrdf_tc [ggrpnr(groups, egcTC, ai)] -= 0.5*imin;
2655 nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5*imin;
2656 if (nrdf_tc[ggrpnr(groups, egcTC, ai)] < 0)
2658 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)]]);
2663 if (ir->nstcomm != 0)
2665 /* Subtract 3 from the number of degrees of freedom in each vcm group
2666 * when com translation is removed and 6 when rotation is removed
2669 switch (ir->comm_mode)
2672 n_sub = ndof_com(ir);
2679 gmx_incons("Checking comm_mode");
2682 for (i = 0; i < groups->grps[egcTC].nr; i++)
2684 /* Count the number of atoms of TC group i for every VCM group */
2685 for (j = 0; j < groups->grps[egcVCM].nr+1; j++)
2690 for (ai = 0; ai < natoms; ai++)
2692 if (ggrpnr(groups, egcTC, ai) == i)
2694 na_vcm[ggrpnr(groups, egcVCM, ai)]++;
2698 /* Correct for VCM removal according to the fraction of each VCM
2699 * group present in this TC group.
2701 nrdf_uc = nrdf_tc[i];
2704 fprintf(debug, "T-group[%d] nrdf_uc = %g, n_sub = %g\n",
2708 for (j = 0; j < groups->grps[egcVCM].nr+1; j++)
2710 if (nrdf_vcm[j] > n_sub)
2712 nrdf_tc[i] += nrdf_uc*((double)na_vcm[j]/(double)na_tot)*
2713 (nrdf_vcm[j] - n_sub)/nrdf_vcm[j];
2717 fprintf(debug, " nrdf_vcm[%d] = %g, nrdf = %g\n",
2718 j, nrdf_vcm[j], nrdf_tc[i]);
2723 for (i = 0; (i < groups->grps[egcTC].nr); i++)
2725 opts->nrdf[i] = nrdf_tc[i];
2726 if (opts->nrdf[i] < 0)
2731 "Number of degrees of freedom in T-Coupling group %s is %.2f\n",
2732 gnames[groups->grps[egcTC].nm_ind[i]], opts->nrdf[i]);
2741 static void decode_cos(char *s, t_cosines *cosine)
2744 char format[STRLEN], f1[STRLEN];
2756 sscanf(t, "%d", &(cosine->n));
2763 snew(cosine->a, cosine->n);
2764 snew(cosine->phi, cosine->n);
2766 sprintf(format, "%%*d");
2767 for (i = 0; (i < cosine->n); i++)
2770 strcat(f1, "%lf%lf");
2771 if (sscanf(t, f1, &a, &phi) < 2)
2773 gmx_fatal(FARGS, "Invalid input for electric field shift: '%s'", t);
2776 cosine->phi[i] = phi;
2777 strcat(format, "%*lf%*lf");
2784 static gmx_bool do_egp_flag(t_inputrec *ir, gmx_groups_t *groups,
2785 const char *option, const char *val, int flag)
2787 /* The maximum number of energy group pairs would be MAXPTR*(MAXPTR+1)/2.
2788 * But since this is much larger than STRLEN, such a line can not be parsed.
2789 * The real maximum is the number of names that fit in a string: STRLEN/2.
2791 #define EGP_MAX (STRLEN/2)
2792 int nelem, i, j, k, nr;
2793 char *names[EGP_MAX];
2797 gnames = groups->grpname;
2799 nelem = str_nelem(val, EGP_MAX, names);
2802 gmx_fatal(FARGS, "The number of groups for %s is odd", option);
2804 nr = groups->grps[egcENER].nr;
2806 for (i = 0; i < nelem/2; i++)
2810 gmx_strcasecmp(names[2*i], *(gnames[groups->grps[egcENER].nm_ind[j]])))
2816 gmx_fatal(FARGS, "%s in %s is not an energy group\n",
2817 names[2*i], option);
2821 gmx_strcasecmp(names[2*i+1], *(gnames[groups->grps[egcENER].nm_ind[k]])))
2827 gmx_fatal(FARGS, "%s in %s is not an energy group\n",
2828 names[2*i+1], option);
2830 if ((j < nr) && (k < nr))
2832 ir->opts.egp_flags[nr*j+k] |= flag;
2833 ir->opts.egp_flags[nr*k+j] |= flag;
2841 void do_index(const char* mdparin, const char *ndx,
2844 t_inputrec *ir, rvec *v,
2848 gmx_groups_t *groups;
2852 char warnbuf[STRLEN], **gnames;
2853 int nr, ntcg, ntau_t, nref_t, nacc, nofg, nSA, nSA_points, nSA_time, nSA_temp;
2856 int nacg, nfreeze, nfrdim, nenergy, nvcm, nuser;
2857 char *ptr1[MAXPTR], *ptr2[MAXPTR], *ptr3[MAXPTR];
2858 int i, j, k, restnm;
2860 gmx_bool bExcl, bTable, bSetTCpar, bAnneal, bRest;
2861 int nQMmethod, nQMbasis, nQMcharge, nQMmult, nbSH, nCASorb, nCASelec,
2862 nSAon, nSAoff, nSAsteps, nQMg, nbOPT, nbTS;
2863 char warn_buf[STRLEN];
2867 fprintf(stderr, "processing index file...\n");
2873 snew(grps->index, 1);
2875 atoms_all = gmx_mtop_global_atoms(mtop);
2876 analyse(&atoms_all, grps, &gnames, FALSE, TRUE);
2877 free_t_atoms(&atoms_all, FALSE);
2881 grps = init_index(ndx, &gnames);
2884 groups = &mtop->groups;
2885 natoms = mtop->natoms;
2886 symtab = &mtop->symtab;
2888 snew(groups->grpname, grps->nr+1);
2890 for (i = 0; (i < grps->nr); i++)
2892 groups->grpname[i] = put_symtab(symtab, gnames[i]);
2894 groups->grpname[i] = put_symtab(symtab, "rest");
2896 srenew(gnames, grps->nr+1);
2897 gnames[restnm] = *(groups->grpname[i]);
2898 groups->ngrpname = grps->nr+1;
2900 set_warning_line(wi, mdparin, -1);
2902 ntau_t = str_nelem(tau_t, MAXPTR, ptr1);
2903 nref_t = str_nelem(ref_t, MAXPTR, ptr2);
2904 ntcg = str_nelem(tcgrps, MAXPTR, ptr3);
2905 if ((ntau_t != ntcg) || (nref_t != ntcg))
2907 gmx_fatal(FARGS, "Invalid T coupling input: %d groups, %d ref-t values and "
2908 "%d tau-t values", ntcg, nref_t, ntau_t);
2911 bSetTCpar = (ir->etc || EI_SD(ir->eI) || ir->eI == eiBD || EI_TPI(ir->eI));
2912 do_numbering(natoms, groups, ntcg, ptr3, grps, gnames, egcTC,
2913 restnm, bSetTCpar ? egrptpALL : egrptpALL_GENREST, bVerbose, wi);
2914 nr = groups->grps[egcTC].nr;
2916 snew(ir->opts.nrdf, nr);
2917 snew(ir->opts.tau_t, nr);
2918 snew(ir->opts.ref_t, nr);
2919 if (ir->eI == eiBD && ir->bd_fric == 0)
2921 fprintf(stderr, "bd-fric=0, so tau-t will be used as the inverse friction constant(s)\n");
2928 gmx_fatal(FARGS, "Not enough ref-t and tau-t values!");
2932 for (i = 0; (i < nr); i++)
2934 ir->opts.tau_t[i] = strtod(ptr1[i], NULL);
2935 if ((ir->eI == eiBD || ir->eI == eiSD2) && ir->opts.tau_t[i] <= 0)
2937 sprintf(warn_buf, "With integrator %s tau-t should be larger than 0", ei_names[ir->eI]);
2938 warning_error(wi, warn_buf);
2941 if (ir->etc != etcVRESCALE && ir->opts.tau_t[i] == 0)
2943 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.");
2946 if (ir->opts.tau_t[i] >= 0)
2948 tau_min = min(tau_min, ir->opts.tau_t[i]);
2951 if (ir->etc != etcNO && ir->nsttcouple == -1)
2953 ir->nsttcouple = ir_optimal_nsttcouple(ir);
2958 if ((ir->etc == etcNOSEHOOVER) && (ir->epc == epcBERENDSEN))
2960 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");
2962 if ((ir->epc == epcMTTK) && (ir->etc > etcNO))
2964 if (ir->nstpcouple != ir->nsttcouple)
2966 int mincouple = min(ir->nstpcouple, ir->nsttcouple);
2967 ir->nstpcouple = ir->nsttcouple = mincouple;
2968 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);
2969 warning_note(wi, warn_buf);
2973 /* velocity verlet with averaged kinetic energy KE = 0.5*(v(t+1/2) - v(t-1/2)) is implemented
2974 primarily for testing purposes, and does not work with temperature coupling other than 1 */
2976 if (ETC_ANDERSEN(ir->etc))
2978 if (ir->nsttcouple != 1)
2981 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");
2982 warning_note(wi, warn_buf);
2985 nstcmin = tcouple_min_integration_steps(ir->etc);
2988 if (tau_min/(ir->delta_t*ir->nsttcouple) < nstcmin)
2990 sprintf(warn_buf, "For proper integration of the %s thermostat, tau-t (%g) should be at least %d times larger than nsttcouple*dt (%g)",
2991 ETCOUPLTYPE(ir->etc),
2993 ir->nsttcouple*ir->delta_t);
2994 warning(wi, warn_buf);
2997 for (i = 0; (i < nr); i++)
2999 ir->opts.ref_t[i] = strtod(ptr2[i], NULL);
3000 if (ir->opts.ref_t[i] < 0)
3002 gmx_fatal(FARGS, "ref-t for group %d negative", i);
3005 /* set the lambda mc temperature to the md integrator temperature (which should be defined
3006 if we are in this conditional) if mc_temp is negative */
3007 if (ir->expandedvals->mc_temp < 0)
3009 ir->expandedvals->mc_temp = ir->opts.ref_t[0]; /*for now, set to the first reft */
3013 /* Simulated annealing for each group. There are nr groups */
3014 nSA = str_nelem(anneal, MAXPTR, ptr1);
3015 if (nSA == 1 && (ptr1[0][0] == 'n' || ptr1[0][0] == 'N'))
3019 if (nSA > 0 && nSA != nr)
3021 gmx_fatal(FARGS, "Not enough annealing values: %d (for %d groups)\n", nSA, nr);
3025 snew(ir->opts.annealing, nr);
3026 snew(ir->opts.anneal_npoints, nr);
3027 snew(ir->opts.anneal_time, nr);
3028 snew(ir->opts.anneal_temp, nr);
3029 for (i = 0; i < nr; i++)
3031 ir->opts.annealing[i] = eannNO;
3032 ir->opts.anneal_npoints[i] = 0;
3033 ir->opts.anneal_time[i] = NULL;
3034 ir->opts.anneal_temp[i] = NULL;
3039 for (i = 0; i < nr; i++)
3041 if (ptr1[i][0] == 'n' || ptr1[i][0] == 'N')
3043 ir->opts.annealing[i] = eannNO;
3045 else if (ptr1[i][0] == 's' || ptr1[i][0] == 'S')
3047 ir->opts.annealing[i] = eannSINGLE;
3050 else if (ptr1[i][0] == 'p' || ptr1[i][0] == 'P')
3052 ir->opts.annealing[i] = eannPERIODIC;
3058 /* Read the other fields too */
3059 nSA_points = str_nelem(anneal_npoints, MAXPTR, ptr1);
3060 if (nSA_points != nSA)
3062 gmx_fatal(FARGS, "Found %d annealing-npoints values for %d groups\n", nSA_points, nSA);
3064 for (k = 0, i = 0; i < nr; i++)
3066 ir->opts.anneal_npoints[i] = strtol(ptr1[i], NULL, 10);
3067 if (ir->opts.anneal_npoints[i] == 1)
3069 gmx_fatal(FARGS, "Please specify at least a start and an end point for annealing\n");
3071 snew(ir->opts.anneal_time[i], ir->opts.anneal_npoints[i]);
3072 snew(ir->opts.anneal_temp[i], ir->opts.anneal_npoints[i]);
3073 k += ir->opts.anneal_npoints[i];
3076 nSA_time = str_nelem(anneal_time, MAXPTR, ptr1);
3079 gmx_fatal(FARGS, "Found %d annealing-time values, wanter %d\n", nSA_time, k);
3081 nSA_temp = str_nelem(anneal_temp, MAXPTR, ptr2);
3084 gmx_fatal(FARGS, "Found %d annealing-temp values, wanted %d\n", nSA_temp, k);
3087 for (i = 0, k = 0; i < nr; i++)
3090 for (j = 0; j < ir->opts.anneal_npoints[i]; j++)
3092 ir->opts.anneal_time[i][j] = strtod(ptr1[k], NULL);
3093 ir->opts.anneal_temp[i][j] = strtod(ptr2[k], NULL);
3096 if (ir->opts.anneal_time[i][0] > (ir->init_t+GMX_REAL_EPS))
3098 gmx_fatal(FARGS, "First time point for annealing > init_t.\n");
3104 if (ir->opts.anneal_time[i][j] < ir->opts.anneal_time[i][j-1])
3106 gmx_fatal(FARGS, "Annealing timepoints out of order: t=%f comes after t=%f\n",
3107 ir->opts.anneal_time[i][j], ir->opts.anneal_time[i][j-1]);
3110 if (ir->opts.anneal_temp[i][j] < 0)
3112 gmx_fatal(FARGS, "Found negative temperature in annealing: %f\n", ir->opts.anneal_temp[i][j]);
3117 /* Print out some summary information, to make sure we got it right */
3118 for (i = 0, k = 0; i < nr; i++)
3120 if (ir->opts.annealing[i] != eannNO)
3122 j = groups->grps[egcTC].nm_ind[i];
3123 fprintf(stderr, "Simulated annealing for group %s: %s, %d timepoints\n",
3124 *(groups->grpname[j]), eann_names[ir->opts.annealing[i]],
3125 ir->opts.anneal_npoints[i]);
3126 fprintf(stderr, "Time (ps) Temperature (K)\n");
3127 /* All terms except the last one */
3128 for (j = 0; j < (ir->opts.anneal_npoints[i]-1); j++)
3130 fprintf(stderr, "%9.1f %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3133 /* Finally the last one */
3134 j = ir->opts.anneal_npoints[i]-1;
3135 if (ir->opts.annealing[i] == eannSINGLE)
3137 fprintf(stderr, "%9.1f- %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3141 fprintf(stderr, "%9.1f %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3142 if (fabs(ir->opts.anneal_temp[i][j]-ir->opts.anneal_temp[i][0]) > GMX_REAL_EPS)
3144 warning_note(wi, "There is a temperature jump when your annealing loops back.\n");
3153 if (ir->ePull != epullNO)
3155 make_pull_groups(ir->pull, pull_grp, grps, gnames);
3160 make_rotation_groups(ir->rot, rot_grp, grps, gnames);
3163 nacc = str_nelem(acc, MAXPTR, ptr1);
3164 nacg = str_nelem(accgrps, MAXPTR, ptr2);
3165 if (nacg*DIM != nacc)
3167 gmx_fatal(FARGS, "Invalid Acceleration input: %d groups and %d acc. values",
3170 do_numbering(natoms, groups, nacg, ptr2, grps, gnames, egcACC,
3171 restnm, egrptpALL_GENREST, bVerbose, wi);
3172 nr = groups->grps[egcACC].nr;
3173 snew(ir->opts.acc, nr);
3174 ir->opts.ngacc = nr;
3176 for (i = k = 0; (i < nacg); i++)
3178 for (j = 0; (j < DIM); j++, k++)
3180 ir->opts.acc[i][j] = strtod(ptr1[k], NULL);
3183 for (; (i < nr); i++)
3185 for (j = 0; (j < DIM); j++)
3187 ir->opts.acc[i][j] = 0;
3191 nfrdim = str_nelem(frdim, MAXPTR, ptr1);
3192 nfreeze = str_nelem(freeze, MAXPTR, ptr2);
3193 if (nfrdim != DIM*nfreeze)
3195 gmx_fatal(FARGS, "Invalid Freezing input: %d groups and %d freeze values",
3198 do_numbering(natoms, groups, nfreeze, ptr2, grps, gnames, egcFREEZE,
3199 restnm, egrptpALL_GENREST, bVerbose, wi);
3200 nr = groups->grps[egcFREEZE].nr;
3201 ir->opts.ngfrz = nr;
3202 snew(ir->opts.nFreeze, nr);
3203 for (i = k = 0; (i < nfreeze); i++)
3205 for (j = 0; (j < DIM); j++, k++)
3207 ir->opts.nFreeze[i][j] = (gmx_strncasecmp(ptr1[k], "Y", 1) == 0);
3208 if (!ir->opts.nFreeze[i][j])
3210 if (gmx_strncasecmp(ptr1[k], "N", 1) != 0)
3212 sprintf(warnbuf, "Please use Y(ES) or N(O) for freezedim only "
3213 "(not %s)", ptr1[k]);
3214 warning(wi, warn_buf);
3219 for (; (i < nr); i++)
3221 for (j = 0; (j < DIM); j++)
3223 ir->opts.nFreeze[i][j] = 0;
3227 nenergy = str_nelem(energy, MAXPTR, ptr1);
3228 do_numbering(natoms, groups, nenergy, ptr1, grps, gnames, egcENER,
3229 restnm, egrptpALL_GENREST, bVerbose, wi);
3230 add_wall_energrps(groups, ir->nwall, symtab);
3231 ir->opts.ngener = groups->grps[egcENER].nr;
3232 nvcm = str_nelem(vcm, MAXPTR, ptr1);
3234 do_numbering(natoms, groups, nvcm, ptr1, grps, gnames, egcVCM,
3235 restnm, nvcm == 0 ? egrptpALL_GENREST : egrptpPART, bVerbose, wi);
3238 warning(wi, "Some atoms are not part of any center of mass motion removal group.\n"
3239 "This may lead to artifacts.\n"
3240 "In most cases one should use one group for the whole system.");
3243 /* Now we have filled the freeze struct, so we can calculate NRDF */
3244 calc_nrdf(mtop, ir, gnames);
3250 /* Must check per group! */
3251 for (i = 0; (i < ir->opts.ngtc); i++)
3253 ntot += ir->opts.nrdf[i];
3255 if (ntot != (DIM*natoms))
3257 fac = sqrt(ntot/(DIM*natoms));
3260 fprintf(stderr, "Scaling velocities by a factor of %.3f to account for constraints\n"
3261 "and removal of center of mass motion\n", fac);
3263 for (i = 0; (i < natoms); i++)
3265 svmul(fac, v[i], v[i]);
3270 nuser = str_nelem(user1, MAXPTR, ptr1);
3271 do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcUser1,
3272 restnm, egrptpALL_GENREST, bVerbose, wi);
3273 nuser = str_nelem(user2, MAXPTR, ptr1);
3274 do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcUser2,
3275 restnm, egrptpALL_GENREST, bVerbose, wi);
3276 nuser = str_nelem(xtc_grps, MAXPTR, ptr1);
3277 do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcXTC,
3278 restnm, egrptpONE, bVerbose, wi);
3279 nofg = str_nelem(orirefitgrp, MAXPTR, ptr1);
3280 do_numbering(natoms, groups, nofg, ptr1, grps, gnames, egcORFIT,
3281 restnm, egrptpALL_GENREST, bVerbose, wi);
3283 /* QMMM input processing */
3284 nQMg = str_nelem(QMMM, MAXPTR, ptr1);
3285 nQMmethod = str_nelem(QMmethod, MAXPTR, ptr2);
3286 nQMbasis = str_nelem(QMbasis, MAXPTR, ptr3);
3287 if ((nQMmethod != nQMg) || (nQMbasis != nQMg))
3289 gmx_fatal(FARGS, "Invalid QMMM input: %d groups %d basissets"
3290 " and %d methods\n", nQMg, nQMbasis, nQMmethod);
3292 /* group rest, if any, is always MM! */
3293 do_numbering(natoms, groups, nQMg, ptr1, grps, gnames, egcQMMM,
3294 restnm, egrptpALL_GENREST, bVerbose, wi);
3295 nr = nQMg; /*atoms->grps[egcQMMM].nr;*/
3296 ir->opts.ngQM = nQMg;
3297 snew(ir->opts.QMmethod, nr);
3298 snew(ir->opts.QMbasis, nr);
3299 for (i = 0; i < nr; i++)
3301 /* input consists of strings: RHF CASSCF PM3 .. These need to be
3302 * converted to the corresponding enum in names.c
3304 ir->opts.QMmethod[i] = search_QMstring(ptr2[i], eQMmethodNR,
3306 ir->opts.QMbasis[i] = search_QMstring(ptr3[i], eQMbasisNR,
3310 nQMmult = str_nelem(QMmult, MAXPTR, ptr1);
3311 nQMcharge = str_nelem(QMcharge, MAXPTR, ptr2);
3312 nbSH = str_nelem(bSH, MAXPTR, ptr3);
3313 snew(ir->opts.QMmult, nr);
3314 snew(ir->opts.QMcharge, nr);
3315 snew(ir->opts.bSH, nr);
3317 for (i = 0; i < nr; i++)
3319 ir->opts.QMmult[i] = strtol(ptr1[i], NULL, 10);
3320 ir->opts.QMcharge[i] = strtol(ptr2[i], NULL, 10);
3321 ir->opts.bSH[i] = (gmx_strncasecmp(ptr3[i], "Y", 1) == 0);
3324 nCASelec = str_nelem(CASelectrons, MAXPTR, ptr1);
3325 nCASorb = str_nelem(CASorbitals, MAXPTR, ptr2);
3326 snew(ir->opts.CASelectrons, nr);
3327 snew(ir->opts.CASorbitals, nr);
3328 for (i = 0; i < nr; i++)
3330 ir->opts.CASelectrons[i] = strtol(ptr1[i], NULL, 10);
3331 ir->opts.CASorbitals[i] = strtol(ptr2[i], NULL, 10);
3333 /* special optimization options */
3335 nbOPT = str_nelem(bOPT, MAXPTR, ptr1);
3336 nbTS = str_nelem(bTS, MAXPTR, ptr2);
3337 snew(ir->opts.bOPT, nr);
3338 snew(ir->opts.bTS, nr);
3339 for (i = 0; i < nr; i++)
3341 ir->opts.bOPT[i] = (gmx_strncasecmp(ptr1[i], "Y", 1) == 0);
3342 ir->opts.bTS[i] = (gmx_strncasecmp(ptr2[i], "Y", 1) == 0);
3344 nSAon = str_nelem(SAon, MAXPTR, ptr1);
3345 nSAoff = str_nelem(SAoff, MAXPTR, ptr2);
3346 nSAsteps = str_nelem(SAsteps, MAXPTR, ptr3);
3347 snew(ir->opts.SAon, nr);
3348 snew(ir->opts.SAoff, nr);
3349 snew(ir->opts.SAsteps, nr);
3351 for (i = 0; i < nr; i++)
3353 ir->opts.SAon[i] = strtod(ptr1[i], NULL);
3354 ir->opts.SAoff[i] = strtod(ptr2[i], NULL);
3355 ir->opts.SAsteps[i] = strtol(ptr3[i], NULL, 10);
3357 /* end of QMMM input */
3361 for (i = 0; (i < egcNR); i++)
3363 fprintf(stderr, "%-16s has %d element(s):", gtypes[i], groups->grps[i].nr);
3364 for (j = 0; (j < groups->grps[i].nr); j++)
3366 fprintf(stderr, " %s", *(groups->grpname[groups->grps[i].nm_ind[j]]));
3368 fprintf(stderr, "\n");
3372 nr = groups->grps[egcENER].nr;
3373 snew(ir->opts.egp_flags, nr*nr);
3375 bExcl = do_egp_flag(ir, groups, "energygrp-excl", egpexcl, EGP_EXCL);
3376 if (bExcl && ir->cutoff_scheme == ecutsVERLET)
3378 warning_error(wi, "Energy group exclusions are not (yet) implemented for the Verlet scheme");
3380 if (bExcl && EEL_FULL(ir->coulombtype))
3382 warning(wi, "Can not exclude the lattice Coulomb energy between energy groups");
3385 bTable = do_egp_flag(ir, groups, "energygrp-table", egptable, EGP_TABLE);
3386 if (bTable && !(ir->vdwtype == evdwUSER) &&
3387 !(ir->coulombtype == eelUSER) && !(ir->coulombtype == eelPMEUSER) &&
3388 !(ir->coulombtype == eelPMEUSERSWITCH))
3390 gmx_fatal(FARGS, "Can only have energy group pair tables in combination with user tables for VdW and/or Coulomb");
3393 decode_cos(efield_x, &(ir->ex[XX]));
3394 decode_cos(efield_xt, &(ir->et[XX]));
3395 decode_cos(efield_y, &(ir->ex[YY]));
3396 decode_cos(efield_yt, &(ir->et[YY]));
3397 decode_cos(efield_z, &(ir->ex[ZZ]));
3398 decode_cos(efield_zt, &(ir->et[ZZ]));
3402 do_adress_index(ir->adress, groups, gnames, &(ir->opts), wi);
3405 for (i = 0; (i < grps->nr); i++)
3417 static void check_disre(gmx_mtop_t *mtop)
3419 gmx_ffparams_t *ffparams;
3420 t_functype *functype;
3422 int i, ndouble, ftype;
3423 int label, old_label;
3425 if (gmx_mtop_ftype_count(mtop, F_DISRES) > 0)
3427 ffparams = &mtop->ffparams;
3428 functype = ffparams->functype;
3429 ip = ffparams->iparams;
3432 for (i = 0; i < ffparams->ntypes; i++)
3434 ftype = functype[i];
3435 if (ftype == F_DISRES)
3437 label = ip[i].disres.label;
3438 if (label == old_label)
3440 fprintf(stderr, "Distance restraint index %d occurs twice\n", label);
3448 gmx_fatal(FARGS, "Found %d double distance restraint indices,\n"
3449 "probably the parameters for multiple pairs in one restraint "
3450 "are not identical\n", ndouble);
3455 static gmx_bool absolute_reference(t_inputrec *ir, gmx_mtop_t *sys,
3456 gmx_bool posres_only,
3460 gmx_mtop_ilistloop_t iloop;
3470 for (d = 0; d < DIM; d++)
3472 AbsRef[d] = (d < ndof_com(ir) ? 0 : 1);
3474 /* Check for freeze groups */
3475 for (g = 0; g < ir->opts.ngfrz; g++)
3477 for (d = 0; d < DIM; d++)
3479 if (ir->opts.nFreeze[g][d] != 0)
3487 /* Check for position restraints */
3488 iloop = gmx_mtop_ilistloop_init(sys);
3489 while (gmx_mtop_ilistloop_next(iloop, &ilist, &nmol))
3492 (AbsRef[XX] == 0 || AbsRef[YY] == 0 || AbsRef[ZZ] == 0))
3494 for (i = 0; i < ilist[F_POSRES].nr; i += 2)
3496 pr = &sys->ffparams.iparams[ilist[F_POSRES].iatoms[i]];
3497 for (d = 0; d < DIM; d++)
3499 if (pr->posres.fcA[d] != 0)
3505 for (i = 0; i < ilist[F_FBPOSRES].nr; i += 2)
3507 /* Check for flat-bottom posres */
3508 pr = &sys->ffparams.iparams[ilist[F_FBPOSRES].iatoms[i]];
3509 if (pr->fbposres.k != 0)
3511 switch (pr->fbposres.geom)
3513 case efbposresSPHERE:
3514 AbsRef[XX] = AbsRef[YY] = AbsRef[ZZ] = 1;
3516 case efbposresCYLINDER:
3517 AbsRef[XX] = AbsRef[YY] = 1;
3519 case efbposresX: /* d=XX */
3520 case efbposresY: /* d=YY */
3521 case efbposresZ: /* d=ZZ */
3522 d = pr->fbposres.geom - efbposresX;
3526 gmx_fatal(FARGS, " Invalid geometry for flat-bottom position restraint.\n"
3527 "Expected nr between 1 and %d. Found %d\n", efbposresNR-1,
3535 return (AbsRef[XX] != 0 && AbsRef[YY] != 0 && AbsRef[ZZ] != 0);
3538 void triple_check(const char *mdparin, t_inputrec *ir, gmx_mtop_t *sys,
3542 int i, m, g, nmol, npct;
3543 gmx_bool bCharge, bAcc;
3544 real gdt_max, *mgrp, mt;
3546 gmx_mtop_atomloop_block_t aloopb;
3547 gmx_mtop_atomloop_all_t aloop;
3550 char warn_buf[STRLEN];
3552 set_warning_line(wi, mdparin, -1);
3554 if (EI_DYNAMICS(ir->eI) && !EI_SD(ir->eI) && ir->eI != eiBD &&
3555 ir->comm_mode == ecmNO &&
3556 !(absolute_reference(ir, sys, FALSE, AbsRef) || ir->nsteps <= 10))
3558 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");
3561 /* Check for pressure coupling with absolute position restraints */
3562 if (ir->epc != epcNO && ir->refcoord_scaling == erscNO)
3564 absolute_reference(ir, sys, TRUE, AbsRef);
3566 for (m = 0; m < DIM; m++)
3568 if (AbsRef[m] && norm2(ir->compress[m]) > 0)
3570 warning(wi, "You are using pressure coupling with absolute position restraints, this will give artifacts. Use the refcoord_scaling option.");
3578 aloopb = gmx_mtop_atomloop_block_init(sys);
3579 while (gmx_mtop_atomloop_block_next(aloopb, &atom, &nmol))
3581 if (atom->q != 0 || atom->qB != 0)
3589 if (EEL_FULL(ir->coulombtype))
3592 "You are using full electrostatics treatment %s for a system without charges.\n"
3593 "This costs a lot of performance for just processing zeros, consider using %s instead.\n",
3594 EELTYPE(ir->coulombtype), EELTYPE(eelCUT));
3595 warning(wi, err_buf);
3600 if (ir->coulombtype == eelCUT && ir->rcoulomb > 0 && !ir->implicit_solvent)
3603 "You are using a plain Coulomb cut-off, which might produce artifacts.\n"
3604 "You might want to consider using %s electrostatics.\n",
3606 warning_note(wi, err_buf);
3610 /* Generalized reaction field */
3611 if (ir->opts.ngtc == 0)
3613 sprintf(err_buf, "No temperature coupling while using coulombtype %s",
3615 CHECK(ir->coulombtype == eelGRF);
3619 sprintf(err_buf, "When using coulombtype = %s"
3620 " ref-t for temperature coupling should be > 0",
3622 CHECK((ir->coulombtype == eelGRF) && (ir->opts.ref_t[0] <= 0));
3625 if (ir->eI == eiSD1 &&
3626 (gmx_mtop_ftype_count(sys, F_CONSTR) > 0 ||
3627 gmx_mtop_ftype_count(sys, F_SETTLE) > 0))
3629 sprintf(warn_buf, "With constraints integrator %s is less accurate, consider using %s instead", ei_names[ir->eI], ei_names[eiSD2]);
3630 warning_note(wi, warn_buf);
3634 for (i = 0; (i < sys->groups.grps[egcACC].nr); i++)
3636 for (m = 0; (m < DIM); m++)
3638 if (fabs(ir->opts.acc[i][m]) > 1e-6)
3647 snew(mgrp, sys->groups.grps[egcACC].nr);
3648 aloop = gmx_mtop_atomloop_all_init(sys);
3649 while (gmx_mtop_atomloop_all_next(aloop, &i, &atom))
3651 mgrp[ggrpnr(&sys->groups, egcACC, i)] += atom->m;
3654 for (i = 0; (i < sys->groups.grps[egcACC].nr); i++)
3656 for (m = 0; (m < DIM); m++)
3658 acc[m] += ir->opts.acc[i][m]*mgrp[i];
3662 for (m = 0; (m < DIM); m++)
3664 if (fabs(acc[m]) > 1e-6)
3666 const char *dim[DIM] = { "X", "Y", "Z" };
3668 "Net Acceleration in %s direction, will %s be corrected\n",
3669 dim[m], ir->nstcomm != 0 ? "" : "not");
3670 if (ir->nstcomm != 0 && m < ndof_com(ir))
3673 for (i = 0; (i < sys->groups.grps[egcACC].nr); i++)
3675 ir->opts.acc[i][m] -= acc[m];
3683 if (ir->efep != efepNO && ir->fepvals->sc_alpha != 0 &&
3684 !gmx_within_tol(sys->ffparams.reppow, 12.0, 10*GMX_DOUBLE_EPS))
3686 gmx_fatal(FARGS, "Soft-core interactions are only supported with VdW repulsion power 12");
3689 if (ir->ePull != epullNO)
3691 if (ir->pull->grp[0].nat == 0)
3693 absolute_reference(ir, sys, FALSE, AbsRef);
3694 for (m = 0; m < DIM; m++)
3696 if (ir->pull->dim[m] && !AbsRef[m])
3698 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.");
3704 if (ir->pull->eGeom == epullgDIRPBC)
3706 for (i = 0; i < 3; i++)
3708 for (m = 0; m <= i; m++)
3710 if ((ir->epc != epcNO && ir->compress[i][m] != 0) ||
3711 ir->deform[i][m] != 0)
3713 for (g = 1; g < ir->pull->ngrp; g++)
3715 if (ir->pull->grp[g].vec[m] != 0)
3717 gmx_fatal(FARGS, "Can not have dynamic box while using pull geometry '%s' (dim %c)", EPULLGEOM(ir->pull->eGeom), 'x'+m);
3729 void double_check(t_inputrec *ir, matrix box, gmx_bool bConstr, warninp_t wi)
3733 char warn_buf[STRLEN];
3736 ptr = check_box(ir->ePBC, box);
3739 warning_error(wi, ptr);
3742 if (bConstr && ir->eConstrAlg == econtSHAKE)
3744 if (ir->shake_tol <= 0.0)
3746 sprintf(warn_buf, "ERROR: shake-tol must be > 0 instead of %g\n",
3748 warning_error(wi, warn_buf);
3751 if (IR_TWINRANGE(*ir) && ir->nstlist > 1)
3753 sprintf(warn_buf, "With twin-range cut-off's and SHAKE the virial and the pressure are incorrect.");
3754 if (ir->epc == epcNO)
3756 warning(wi, warn_buf);
3760 warning_error(wi, warn_buf);
3765 if ( (ir->eConstrAlg == econtLINCS) && bConstr)
3767 /* If we have Lincs constraints: */
3768 if (ir->eI == eiMD && ir->etc == etcNO &&
3769 ir->eConstrAlg == econtLINCS && ir->nLincsIter == 1)
3771 sprintf(warn_buf, "For energy conservation with LINCS, lincs_iter should be 2 or larger.\n");
3772 warning_note(wi, warn_buf);
3775 if ((ir->eI == eiCG || ir->eI == eiLBFGS) && (ir->nProjOrder < 8))
3777 sprintf(warn_buf, "For accurate %s with LINCS constraints, lincs-order should be 8 or more.", ei_names[ir->eI]);
3778 warning_note(wi, warn_buf);
3780 if (ir->epc == epcMTTK)
3782 warning_error(wi, "MTTK not compatible with lincs -- use shake instead.");
3786 if (ir->LincsWarnAngle > 90.0)
3788 sprintf(warn_buf, "lincs-warnangle can not be larger than 90 degrees, setting it to 90.\n");
3789 warning(wi, warn_buf);
3790 ir->LincsWarnAngle = 90.0;
3793 if (ir->ePBC != epbcNONE)
3795 if (ir->nstlist == 0)
3797 warning(wi, "With nstlist=0 atoms are only put into the box at step 0, therefore drifting atoms might cause the simulation to crash.");
3799 bTWIN = (ir->rlistlong > ir->rlist);
3800 if (ir->ns_type == ensGRID)
3802 if (sqr(ir->rlistlong) >= max_cutoff2(ir->ePBC, box))
3804 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",
3805 bTWIN ? (ir->rcoulomb == ir->rlistlong ? "rcoulomb" : "rvdw") : "rlist");
3806 warning_error(wi, warn_buf);
3811 min_size = min(box[XX][XX], min(box[YY][YY], box[ZZ][ZZ]));
3812 if (2*ir->rlistlong >= min_size)
3814 sprintf(warn_buf, "ERROR: One of the box lengths is smaller than twice the cut-off length. Increase the box size or decrease rlist.");
3815 warning_error(wi, warn_buf);
3818 fprintf(stderr, "Grid search might allow larger cut-off's than simple search with triclinic boxes.");
3825 void check_chargegroup_radii(const gmx_mtop_t *mtop, const t_inputrec *ir,
3829 real rvdw1, rvdw2, rcoul1, rcoul2;
3830 char warn_buf[STRLEN];
3832 calc_chargegroup_radii(mtop, x, &rvdw1, &rvdw2, &rcoul1, &rcoul2);
3836 printf("Largest charge group radii for Van der Waals: %5.3f, %5.3f nm\n",
3841 printf("Largest charge group radii for Coulomb: %5.3f, %5.3f nm\n",
3847 if (rvdw1 + rvdw2 > ir->rlist ||
3848 rcoul1 + rcoul2 > ir->rlist)
3851 "The sum of the two largest charge group radii (%f) "
3852 "is larger than rlist (%f)\n",
3853 max(rvdw1+rvdw2, rcoul1+rcoul2), ir->rlist);
3854 warning(wi, warn_buf);
3858 /* Here we do not use the zero at cut-off macro,
3859 * since user defined interactions might purposely
3860 * not be zero at the cut-off.
3862 if ((EVDW_IS_ZERO_AT_CUTOFF(ir->vdwtype) ||
3863 ir->vdw_modifier != eintmodNONE) &&
3864 rvdw1 + rvdw2 > ir->rlistlong - ir->rvdw)
3866 sprintf(warn_buf, "The sum of the two largest charge group "
3867 "radii (%f) is larger than %s (%f) - rvdw (%f).\n"
3868 "With exact cut-offs, better performance can be "
3869 "obtained with cutoff-scheme = %s, because it "
3870 "does not use charge groups at all.",
3872 ir->rlistlong > ir->rlist ? "rlistlong" : "rlist",
3873 ir->rlistlong, ir->rvdw,
3874 ecutscheme_names[ecutsVERLET]);
3877 warning(wi, warn_buf);
3881 warning_note(wi, warn_buf);
3884 if ((EEL_IS_ZERO_AT_CUTOFF(ir->coulombtype) ||
3885 ir->coulomb_modifier != eintmodNONE) &&
3886 rcoul1 + rcoul2 > ir->rlistlong - ir->rcoulomb)
3888 sprintf(warn_buf, "The sum of the two largest charge group radii (%f) is larger than %s (%f) - rcoulomb (%f).\n"
3889 "With exact cut-offs, better performance can be obtained with cutoff-scheme = %s, because it does not use charge groups at all.",
3891 ir->rlistlong > ir->rlist ? "rlistlong" : "rlist",
3892 ir->rlistlong, ir->rcoulomb,
3893 ecutscheme_names[ecutsVERLET]);
3896 warning(wi, warn_buf);
3900 warning_note(wi, warn_buf);