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. */
104 void init_ir(t_inputrec *ir, t_gromppopts *opts)
106 snew(opts->include, STRLEN);
107 snew(opts->define, STRLEN);
108 snew(ir->fepvals, 1);
109 snew(ir->expandedvals, 1);
110 snew(ir->simtempvals, 1);
113 static void GetSimTemps(int ntemps, t_simtemp *simtemp, double *temperature_lambdas)
118 for (i = 0; i < ntemps; i++)
120 /* simple linear scaling -- allows more control */
121 if (simtemp->eSimTempScale == esimtempLINEAR)
123 simtemp->temperatures[i] = simtemp->simtemp_low + (simtemp->simtemp_high-simtemp->simtemp_low)*temperature_lambdas[i];
125 else if (simtemp->eSimTempScale == esimtempGEOMETRIC) /* should give roughly equal acceptance for constant heat capacity . . . */
127 simtemp->temperatures[i] = simtemp->simtemp_low * pow(simtemp->simtemp_high/simtemp->simtemp_low, (1.0*i)/(ntemps-1));
129 else if (simtemp->eSimTempScale == esimtempEXPONENTIAL)
131 simtemp->temperatures[i] = simtemp->simtemp_low + (simtemp->simtemp_high-simtemp->simtemp_low)*((exp(temperature_lambdas[i])-1)/(exp(1.0)-1));
136 sprintf(errorstr, "eSimTempScale=%d not defined", simtemp->eSimTempScale);
137 gmx_fatal(FARGS, errorstr);
144 static void _low_check(gmx_bool b, char *s, warninp_t wi)
148 warning_error(wi, s);
152 static void check_nst(const char *desc_nst, int nst,
153 const char *desc_p, int *p,
158 if (*p > 0 && *p % nst != 0)
160 /* Round up to the next multiple of nst */
161 *p = ((*p)/nst + 1)*nst;
162 sprintf(buf, "%s should be a multiple of %s, changing %s to %d\n",
163 desc_p, desc_nst, desc_p, *p);
168 static gmx_bool ir_NVE(const t_inputrec *ir)
170 return ((ir->eI == eiMD || EI_VV(ir->eI)) && ir->etc == etcNO);
173 static int lcd(int n1, int n2)
178 for (i = 2; (i <= n1 && i <= n2); i++)
180 if (n1 % i == 0 && n2 % i == 0)
189 static void process_interaction_modifier(const t_inputrec *ir, int *eintmod)
191 if (*eintmod == eintmodPOTSHIFT_VERLET)
193 if (ir->cutoff_scheme == ecutsVERLET)
195 *eintmod = eintmodPOTSHIFT;
199 *eintmod = eintmodNONE;
204 void check_ir(const char *mdparin, t_inputrec *ir, t_gromppopts *opts,
206 /* Check internal consistency */
208 /* Strange macro: first one fills the err_buf, and then one can check
209 * the condition, which will print the message and increase the error
212 #define CHECK(b) _low_check(b, err_buf, wi)
213 char err_buf[256], warn_buf[STRLEN];
219 t_lambda *fep = ir->fepvals;
220 t_expanded *expand = ir->expandedvals;
222 set_warning_line(wi, mdparin, -1);
224 /* BASIC CUT-OFF STUFF */
225 if (ir->rcoulomb < 0)
227 warning_error(wi, "rcoulomb should be >= 0");
231 warning_error(wi, "rvdw should be >= 0");
234 !(ir->cutoff_scheme == ecutsVERLET && ir->verletbuf_drift > 0))
236 warning_error(wi, "rlist should be >= 0");
239 process_interaction_modifier(ir, &ir->coulomb_modifier);
240 process_interaction_modifier(ir, &ir->vdw_modifier);
242 if (ir->cutoff_scheme == ecutsGROUP)
244 /* BASIC CUT-OFF STUFF */
245 if (ir->rlist == 0 ||
246 !((EEL_MIGHT_BE_ZERO_AT_CUTOFF(ir->coulombtype) && ir->rcoulomb > ir->rlist) ||
247 (EVDW_MIGHT_BE_ZERO_AT_CUTOFF(ir->vdwtype) && ir->rvdw > ir->rlist)))
249 /* No switched potential and/or no twin-range:
250 * we can set the long-range cut-off to the maximum of the other cut-offs.
252 ir->rlistlong = max_cutoff(ir->rlist, max_cutoff(ir->rvdw, ir->rcoulomb));
254 else if (ir->rlistlong < 0)
256 ir->rlistlong = max_cutoff(ir->rlist, max_cutoff(ir->rvdw, ir->rcoulomb));
257 sprintf(warn_buf, "rlistlong was not set, setting it to %g (no buffer)",
259 warning(wi, warn_buf);
261 if (ir->rlistlong == 0 && ir->ePBC != epbcNONE)
263 warning_error(wi, "Can not have an infinite cut-off with PBC");
265 if (ir->rlistlong > 0 && (ir->rlist == 0 || ir->rlistlong < ir->rlist))
267 warning_error(wi, "rlistlong can not be shorter than rlist");
269 if (IR_TWINRANGE(*ir) && ir->nstlist <= 0)
271 warning_error(wi, "Can not have nstlist<=0 with twin-range interactions");
275 if (ir->rlistlong == ir->rlist)
279 else if (ir->rlistlong > ir->rlist && ir->nstcalclr == 0)
281 warning_error(wi, "With different cutoffs for electrostatics and VdW, nstcalclr must be -1 or a positive number");
284 if (ir->cutoff_scheme == ecutsVERLET)
288 /* Normal Verlet type neighbor-list, currently only limited feature support */
289 if (inputrec2nboundeddim(ir) < 3)
291 warning_error(wi, "With Verlet lists only full pbc or pbc=xy with walls is supported");
293 if (ir->rcoulomb != ir->rvdw)
295 warning_error(wi, "With Verlet lists rcoulomb!=rvdw is not supported");
297 if (ir->vdwtype != evdwCUT)
299 warning_error(wi, "With Verlet lists only cut-off LJ interactions are supported");
301 if (!(ir->coulombtype == eelCUT ||
302 (EEL_RF(ir->coulombtype) && ir->coulombtype != eelRF_NEC) ||
303 EEL_PME(ir->coulombtype) || ir->coulombtype == eelEWALD))
305 warning_error(wi, "With Verlet lists only cut-off, reaction-field, PME and Ewald electrostatics are supported");
308 if (ir->nstlist <= 0)
310 warning_error(wi, "With Verlet lists nstlist should be larger than 0");
313 if (ir->nstlist < 10)
315 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.");
318 rc_max = max(ir->rvdw, ir->rcoulomb);
320 if (ir->verletbuf_drift <= 0)
322 if (ir->verletbuf_drift == 0)
324 warning_error(wi, "Can not have an energy drift of exactly 0");
327 if (ir->rlist < rc_max)
329 warning_error(wi, "With verlet lists rlist can not be smaller than rvdw or rcoulomb");
332 if (ir->rlist == rc_max && ir->nstlist > 1)
334 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.");
339 if (ir->rlist > rc_max)
341 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.");
344 if (ir->nstlist == 1)
346 /* No buffer required */
351 if (EI_DYNAMICS(ir->eI))
353 if (EI_MD(ir->eI) && ir->etc == etcNO)
355 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.");
358 if (inputrec2nboundeddim(ir) < 3)
360 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.");
362 /* Set rlist temporarily so we can continue processing */
367 /* Set the buffer to 5% of the cut-off */
368 ir->rlist = 1.05*rc_max;
373 /* No twin-range calculations with Verlet lists */
374 ir->rlistlong = ir->rlist;
377 if (ir->nstcalclr == -1)
379 /* if rlist=rlistlong, this will later be changed to nstcalclr=0 */
380 ir->nstcalclr = ir->nstlist;
382 else if (ir->nstcalclr > 0)
384 if (ir->nstlist > 0 && (ir->nstlist % ir->nstcalclr != 0))
386 warning_error(wi, "nstlist must be evenly divisible by nstcalclr. Use nstcalclr = -1 to automatically follow nstlist");
389 else if (ir->nstcalclr < -1)
391 warning_error(wi, "nstcalclr must be a positive number (divisor of nstcalclr), or -1 to follow nstlist.");
394 if (EEL_PME(ir->coulombtype) && ir->rcoulomb > ir->rvdw && ir->nstcalclr > 1)
396 warning_error(wi, "When used with PME, the long-range component of twin-range interactions must be updated every step (nstcalclr)");
399 /* GENERAL INTEGRATOR STUFF */
400 if (!(ir->eI == eiMD || EI_VV(ir->eI)))
404 if (ir->eI == eiVVAK)
406 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]);
407 warning_note(wi, warn_buf);
409 if (!EI_DYNAMICS(ir->eI))
413 if (EI_DYNAMICS(ir->eI))
415 if (ir->nstcalcenergy < 0)
417 ir->nstcalcenergy = ir_optimal_nstcalcenergy(ir);
418 if (ir->nstenergy != 0 && ir->nstenergy < ir->nstcalcenergy)
420 /* nstcalcenergy larger than nstener does not make sense.
421 * We ideally want nstcalcenergy=nstener.
425 ir->nstcalcenergy = lcd(ir->nstenergy, ir->nstlist);
429 ir->nstcalcenergy = ir->nstenergy;
433 else if ( (ir->nstenergy > 0 && ir->nstcalcenergy > ir->nstenergy) ||
434 (ir->efep != efepNO && ir->fepvals->nstdhdl > 0 &&
435 (ir->nstcalcenergy > ir->fepvals->nstdhdl) ) )
438 const char *nsten = "nstenergy";
439 const char *nstdh = "nstdhdl";
440 const char *min_name = nsten;
441 int min_nst = ir->nstenergy;
443 /* find the smallest of ( nstenergy, nstdhdl ) */
444 if (ir->efep != efepNO && ir->fepvals->nstdhdl > 0 &&
445 (ir->fepvals->nstdhdl < ir->nstenergy) )
447 min_nst = ir->fepvals->nstdhdl;
450 /* If the user sets nstenergy small, we should respect that */
452 "Setting nstcalcenergy (%d) equal to %s (%d)",
453 ir->nstcalcenergy, min_name, min_nst);
454 warning_note(wi, warn_buf);
455 ir->nstcalcenergy = min_nst;
458 if (ir->epc != epcNO)
460 if (ir->nstpcouple < 0)
462 ir->nstpcouple = ir_optimal_nstpcouple(ir);
465 if (IR_TWINRANGE(*ir))
467 check_nst("nstlist", ir->nstlist,
468 "nstcalcenergy", &ir->nstcalcenergy, wi);
469 if (ir->epc != epcNO)
471 check_nst("nstlist", ir->nstlist,
472 "nstpcouple", &ir->nstpcouple, wi);
476 if (ir->nstcalcenergy > 0)
478 if (ir->efep != efepNO)
480 /* nstdhdl should be a multiple of nstcalcenergy */
481 check_nst("nstcalcenergy", ir->nstcalcenergy,
482 "nstdhdl", &ir->fepvals->nstdhdl, wi);
483 /* nstexpanded should be a multiple of nstcalcenergy */
484 check_nst("nstcalcenergy", ir->nstcalcenergy,
485 "nstexpanded", &ir->expandedvals->nstexpanded, wi);
487 /* for storing exact averages nstenergy should be
488 * a multiple of nstcalcenergy
490 check_nst("nstcalcenergy", ir->nstcalcenergy,
491 "nstenergy", &ir->nstenergy, wi);
496 if ((EI_SD(ir->eI) || ir->eI == eiBD) &&
497 ir->bContinuation && ir->ld_seed != -1)
499 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)");
505 sprintf(err_buf, "TPI only works with pbc = %s", epbc_names[epbcXYZ]);
506 CHECK(ir->ePBC != epbcXYZ);
507 sprintf(err_buf, "TPI only works with ns = %s", ens_names[ensGRID]);
508 CHECK(ir->ns_type != ensGRID);
509 sprintf(err_buf, "with TPI nstlist should be larger than zero");
510 CHECK(ir->nstlist <= 0);
511 sprintf(err_buf, "TPI does not work with full electrostatics other than PME");
512 CHECK(EEL_FULL(ir->coulombtype) && !EEL_PME(ir->coulombtype));
516 if ( (opts->nshake > 0) && (opts->bMorse) )
519 "Using morse bond-potentials while constraining bonds is useless");
520 warning(wi, warn_buf);
523 if ((EI_SD(ir->eI) || ir->eI == eiBD) &&
524 ir->bContinuation && ir->ld_seed != -1)
526 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)");
528 /* verify simulated tempering options */
532 gmx_bool bAllTempZero = TRUE;
533 for (i = 0; i < fep->n_lambda; i++)
535 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]);
536 CHECK((fep->all_lambda[efptTEMPERATURE][i] < 0) || (fep->all_lambda[efptTEMPERATURE][i] > 1));
537 if (fep->all_lambda[efptTEMPERATURE][i] > 0)
539 bAllTempZero = FALSE;
542 sprintf(err_buf, "if simulated tempering is on, temperature-lambdas may not be all zero");
543 CHECK(bAllTempZero == TRUE);
545 sprintf(err_buf, "Simulated tempering is currently only compatible with md-vv");
546 CHECK(ir->eI != eiVV);
548 /* check compatability of the temperature coupling with simulated tempering */
550 if (ir->etc == etcNOSEHOOVER)
552 sprintf(warn_buf, "Nose-Hoover based temperature control such as [%s] my not be entirelyconsistent with simulated tempering", etcoupl_names[ir->etc]);
553 warning_note(wi, warn_buf);
556 /* check that the temperatures make sense */
558 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);
559 CHECK(ir->simtempvals->simtemp_high <= ir->simtempvals->simtemp_low);
561 sprintf(err_buf, "Higher simulated tempering temperature (%g) must be >= zero", ir->simtempvals->simtemp_high);
562 CHECK(ir->simtempvals->simtemp_high <= 0);
564 sprintf(err_buf, "Lower simulated tempering temperature (%g) must be >= zero", ir->simtempvals->simtemp_low);
565 CHECK(ir->simtempvals->simtemp_low <= 0);
568 /* verify free energy options */
570 if (ir->efep != efepNO)
573 sprintf(err_buf, "The soft-core power is %d and can only be 1 or 2",
575 CHECK(fep->sc_alpha != 0 && fep->sc_power != 1 && fep->sc_power != 2);
577 sprintf(err_buf, "The soft-core sc-r-power is %d and can only be 6 or 48",
578 (int)fep->sc_r_power);
579 CHECK(fep->sc_alpha != 0 && fep->sc_r_power != 6.0 && fep->sc_r_power != 48.0);
581 sprintf(err_buf, "Can't use postive delta-lambda (%g) if initial state/lambda does not start at zero", fep->delta_lambda);
582 CHECK(fep->delta_lambda > 0 && ((fep->init_fep_state > 0) || (fep->init_lambda > 0)));
584 sprintf(err_buf, "Can't use postive delta-lambda (%g) with expanded ensemble simulations", fep->delta_lambda);
585 CHECK(fep->delta_lambda > 0 && (ir->efep == efepEXPANDED));
587 sprintf(err_buf, "Free-energy not implemented for Ewald");
588 CHECK(ir->coulombtype == eelEWALD);
590 /* check validty of lambda inputs */
591 if (fep->n_lambda == 0)
593 /* Clear output in case of no states:*/
594 sprintf(err_buf, "init-lambda-state set to %d: no lambda states are defined.", fep->init_fep_state);
595 CHECK((fep->init_fep_state >= 0) && (fep->n_lambda == 0));
599 sprintf(err_buf, "initial thermodynamic state %d does not exist, only goes to %d", fep->init_fep_state, fep->n_lambda-1);
600 CHECK((fep->init_fep_state >= fep->n_lambda));
603 sprintf(err_buf, "Lambda state must be set, either with init-lambda-state or with init-lambda");
604 CHECK((fep->init_fep_state < 0) && (fep->init_lambda < 0));
606 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",
607 fep->init_lambda, fep->init_fep_state);
608 CHECK((fep->init_fep_state >= 0) && (fep->init_lambda >= 0));
612 if ((fep->init_lambda >= 0) && (fep->delta_lambda == 0))
616 for (i = 0; i < efptNR; i++)
618 if (fep->separate_dvdl[i])
623 if (n_lambda_terms > 1)
625 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.");
626 warning(wi, warn_buf);
629 if (n_lambda_terms < 2 && fep->n_lambda > 0)
632 "init-lambda is deprecated for setting lambda state (except for slow growth). Use init-lambda-state instead.");
636 for (j = 0; j < efptNR; j++)
638 for (i = 0; i < fep->n_lambda; i++)
640 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]);
641 CHECK((fep->all_lambda[j][i] < 0) || (fep->all_lambda[j][i] > 1));
645 if ((fep->sc_alpha > 0) && (!fep->bScCoul))
647 for (i = 0; i < fep->n_lambda; i++)
649 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],
650 fep->all_lambda[efptCOUL][i]);
651 CHECK((fep->sc_alpha > 0) &&
652 (((fep->all_lambda[efptCOUL][i] > 0.0) &&
653 (fep->all_lambda[efptCOUL][i] < 1.0)) &&
654 ((fep->all_lambda[efptVDW][i] > 0.0) &&
655 (fep->all_lambda[efptVDW][i] < 1.0))));
659 if ((fep->bScCoul) && (EEL_PME(ir->coulombtype)))
661 real sigma, lambda, r_sc;
664 /* Maximum estimate for A and B charges equal with lambda power 1 */
666 r_sc = pow(lambda*fep->sc_alpha*pow(sigma/ir->rcoulomb, fep->sc_r_power) + 1.0, 1.0/fep->sc_r_power);
667 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.",
669 sigma, lambda, r_sc - 1.0, ir->ewald_rtol);
670 warning_note(wi, warn_buf);
673 /* Free Energy Checks -- In an ideal world, slow growth and FEP would
674 be treated differently, but that's the next step */
676 for (i = 0; i < efptNR; i++)
678 for (j = 0; j < fep->n_lambda; j++)
680 sprintf(err_buf, "%s[%d] must be between 0 and 1", efpt_names[i], j);
681 CHECK((fep->all_lambda[i][j] < 0) || (fep->all_lambda[i][j] > 1));
686 if ((ir->bSimTemp) || (ir->efep == efepEXPANDED))
689 expand = ir->expandedvals;
691 /* checking equilibration of weights inputs for validity */
693 sprintf(err_buf, "weight-equil-number-all-lambda (%d) is ignored if lmc-weights-equil is not equal to %s",
694 expand->equil_n_at_lam, elmceq_names[elmceqNUMATLAM]);
695 CHECK((expand->equil_n_at_lam > 0) && (expand->elmceq != elmceqNUMATLAM));
697 sprintf(err_buf, "weight-equil-number-samples (%d) is ignored if lmc-weights-equil is not equal to %s",
698 expand->equil_samples, elmceq_names[elmceqSAMPLES]);
699 CHECK((expand->equil_samples > 0) && (expand->elmceq != elmceqSAMPLES));
701 sprintf(err_buf, "weight-equil-number-steps (%d) is ignored if lmc-weights-equil is not equal to %s",
702 expand->equil_steps, elmceq_names[elmceqSTEPS]);
703 CHECK((expand->equil_steps > 0) && (expand->elmceq != elmceqSTEPS));
705 sprintf(err_buf, "weight-equil-wl-delta (%d) is ignored if lmc-weights-equil is not equal to %s",
706 expand->equil_samples, elmceq_names[elmceqWLDELTA]);
707 CHECK((expand->equil_wl_delta > 0) && (expand->elmceq != elmceqWLDELTA));
709 sprintf(err_buf, "weight-equil-count-ratio (%f) is ignored if lmc-weights-equil is not equal to %s",
710 expand->equil_ratio, elmceq_names[elmceqRATIO]);
711 CHECK((expand->equil_ratio > 0) && (expand->elmceq != elmceqRATIO));
713 sprintf(err_buf, "weight-equil-number-all-lambda (%d) must be a positive integer if lmc-weights-equil=%s",
714 expand->equil_n_at_lam, elmceq_names[elmceqNUMATLAM]);
715 CHECK((expand->equil_n_at_lam <= 0) && (expand->elmceq == elmceqNUMATLAM));
717 sprintf(err_buf, "weight-equil-number-samples (%d) must be a positive integer if lmc-weights-equil=%s",
718 expand->equil_samples, elmceq_names[elmceqSAMPLES]);
719 CHECK((expand->equil_samples <= 0) && (expand->elmceq == elmceqSAMPLES));
721 sprintf(err_buf, "weight-equil-number-steps (%d) must be a positive integer if lmc-weights-equil=%s",
722 expand->equil_steps, elmceq_names[elmceqSTEPS]);
723 CHECK((expand->equil_steps <= 0) && (expand->elmceq == elmceqSTEPS));
725 sprintf(err_buf, "weight-equil-wl-delta (%f) must be > 0 if lmc-weights-equil=%s",
726 expand->equil_wl_delta, elmceq_names[elmceqWLDELTA]);
727 CHECK((expand->equil_wl_delta <= 0) && (expand->elmceq == elmceqWLDELTA));
729 sprintf(err_buf, "weight-equil-count-ratio (%f) must be > 0 if lmc-weights-equil=%s",
730 expand->equil_ratio, elmceq_names[elmceqRATIO]);
731 CHECK((expand->equil_ratio <= 0) && (expand->elmceq == elmceqRATIO));
733 sprintf(err_buf, "lmc-weights-equil=%s only possible when lmc-stats = %s or lmc-stats %s",
734 elmceq_names[elmceqWLDELTA], elamstats_names[elamstatsWL], elamstats_names[elamstatsWWL]);
735 CHECK((expand->elmceq == elmceqWLDELTA) && (!EWL(expand->elamstats)));
737 sprintf(err_buf, "lmc-repeats (%d) must be greater than 0", expand->lmc_repeats);
738 CHECK((expand->lmc_repeats <= 0));
739 sprintf(err_buf, "minimum-var-min (%d) must be greater than 0", expand->minvarmin);
740 CHECK((expand->minvarmin <= 0));
741 sprintf(err_buf, "weight-c-range (%d) must be greater or equal to 0", expand->c_range);
742 CHECK((expand->c_range < 0));
743 sprintf(err_buf, "init-lambda-state (%d) must be zero if lmc-forced-nstart (%d)> 0 and lmc-move != 'no'",
744 fep->init_fep_state, expand->lmc_forced_nstart);
745 CHECK((fep->init_fep_state != 0) && (expand->lmc_forced_nstart > 0) && (expand->elmcmove != elmcmoveNO));
746 sprintf(err_buf, "lmc-forced-nstart (%d) must not be negative", expand->lmc_forced_nstart);
747 CHECK((expand->lmc_forced_nstart < 0));
748 sprintf(err_buf, "init-lambda-state (%d) must be in the interval [0,number of lambdas)", fep->init_fep_state);
749 CHECK((fep->init_fep_state < 0) || (fep->init_fep_state >= fep->n_lambda));
751 sprintf(err_buf, "init-wl-delta (%f) must be greater than or equal to 0", expand->init_wl_delta);
752 CHECK((expand->init_wl_delta < 0));
753 sprintf(err_buf, "wl-ratio (%f) must be between 0 and 1", expand->wl_ratio);
754 CHECK((expand->wl_ratio <= 0) || (expand->wl_ratio >= 1));
755 sprintf(err_buf, "wl-scale (%f) must be between 0 and 1", expand->wl_scale);
756 CHECK((expand->wl_scale <= 0) || (expand->wl_scale >= 1));
758 /* if there is no temperature control, we need to specify an MC temperature */
759 sprintf(err_buf, "If there is no temperature control, and lmc-mcmove!= 'no',mc_temperature must be set to a positive number");
760 if (expand->nstTij > 0)
762 sprintf(err_buf, "nst-transition-matrix (%d) must be an integer multiple of nstlog (%d)",
763 expand->nstTij, ir->nstlog);
764 CHECK((mod(expand->nstTij, ir->nstlog) != 0));
769 sprintf(err_buf, "walls only work with pbc=%s", epbc_names[epbcXY]);
770 CHECK(ir->nwall && ir->ePBC != epbcXY);
773 if (ir->ePBC != epbcXYZ && ir->nwall != 2)
775 if (ir->ePBC == epbcNONE)
777 if (ir->epc != epcNO)
779 warning(wi, "Turning off pressure coupling for vacuum system");
785 sprintf(err_buf, "Can not have pressure coupling with pbc=%s",
786 epbc_names[ir->ePBC]);
787 CHECK(ir->epc != epcNO);
789 sprintf(err_buf, "Can not have Ewald with pbc=%s", epbc_names[ir->ePBC]);
790 CHECK(EEL_FULL(ir->coulombtype));
792 sprintf(err_buf, "Can not have dispersion correction with pbc=%s",
793 epbc_names[ir->ePBC]);
794 CHECK(ir->eDispCorr != edispcNO);
797 if (ir->rlist == 0.0)
799 sprintf(err_buf, "can only have neighborlist cut-off zero (=infinite)\n"
800 "with coulombtype = %s or coulombtype = %s\n"
801 "without periodic boundary conditions (pbc = %s) and\n"
802 "rcoulomb and rvdw set to zero",
803 eel_names[eelCUT], eel_names[eelUSER], epbc_names[epbcNONE]);
804 CHECK(((ir->coulombtype != eelCUT) && (ir->coulombtype != eelUSER)) ||
805 (ir->ePBC != epbcNONE) ||
806 (ir->rcoulomb != 0.0) || (ir->rvdw != 0.0));
810 warning_error(wi, "Can not have heuristic neighborlist updates without cut-off");
814 warning_note(wi, "Simulating without cut-offs is usually (slightly) faster with nstlist=0, nstype=simple and particle decomposition");
819 if (ir->nstcomm == 0)
821 ir->comm_mode = ecmNO;
823 if (ir->comm_mode != ecmNO)
827 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");
828 ir->nstcomm = abs(ir->nstcomm);
831 if (ir->nstcalcenergy > 0 && ir->nstcomm < ir->nstcalcenergy)
833 warning_note(wi, "nstcomm < nstcalcenergy defeats the purpose of nstcalcenergy, setting nstcomm to nstcalcenergy");
834 ir->nstcomm = ir->nstcalcenergy;
837 if (ir->comm_mode == ecmANGULAR)
839 sprintf(err_buf, "Can not remove the rotation around the center of mass with periodic molecules");
840 CHECK(ir->bPeriodicMols);
841 if (ir->ePBC != epbcNONE)
843 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).");
848 if (EI_STATE_VELOCITY(ir->eI) && ir->ePBC == epbcNONE && ir->comm_mode != ecmANGULAR)
850 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.");
853 sprintf(err_buf, "Twin-range neighbour searching (NS) with simple NS"
854 " algorithm not implemented");
855 CHECK(((ir->rcoulomb > ir->rlist) || (ir->rvdw > ir->rlist))
856 && (ir->ns_type == ensSIMPLE));
858 /* TEMPERATURE COUPLING */
859 if (ir->etc == etcYES)
861 ir->etc = etcBERENDSEN;
862 warning_note(wi, "Old option for temperature coupling given: "
863 "changing \"yes\" to \"Berendsen\"\n");
866 if ((ir->etc == etcNOSEHOOVER) || (ir->epc == epcMTTK))
868 if (ir->opts.nhchainlength < 1)
870 sprintf(warn_buf, "number of Nose-Hoover chains (currently %d) cannot be less than 1,reset to 1\n", ir->opts.nhchainlength);
871 ir->opts.nhchainlength = 1;
872 warning(wi, warn_buf);
875 if (ir->etc == etcNOSEHOOVER && !EI_VV(ir->eI) && ir->opts.nhchainlength > 1)
877 warning_note(wi, "leapfrog does not yet support Nose-Hoover chains, nhchainlength reset to 1");
878 ir->opts.nhchainlength = 1;
883 ir->opts.nhchainlength = 0;
886 if (ir->eI == eiVVAK)
888 sprintf(err_buf, "%s implemented primarily for validation, and requires nsttcouple = 1 and nstpcouple = 1.",
890 CHECK((ir->nsttcouple != 1) || (ir->nstpcouple != 1));
893 if (ETC_ANDERSEN(ir->etc))
895 sprintf(err_buf, "%s temperature control not supported for integrator %s.", etcoupl_names[ir->etc], ei_names[ir->eI]);
896 CHECK(!(EI_VV(ir->eI)));
898 for (i = 0; i < ir->opts.ngtc; i++)
900 sprintf(err_buf, "all tau_t must currently be equal using Andersen temperature control, violated for group %d", i);
901 CHECK(ir->opts.tau_t[0] != ir->opts.tau_t[i]);
902 sprintf(err_buf, "all tau_t must be postive using Andersen temperature control, tau_t[%d]=%10.6f",
903 i, ir->opts.tau_t[i]);
904 CHECK(ir->opts.tau_t[i] < 0);
906 if (ir->nstcomm > 0 && (ir->etc == etcANDERSEN))
908 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]);
909 warning_note(wi, warn_buf);
912 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]);
913 CHECK(ir->nstcomm > 1 && (ir->etc == etcANDERSEN));
915 for (i = 0; i < ir->opts.ngtc; i++)
917 int nsteps = (int)(ir->opts.tau_t[i]/ir->delta_t);
918 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);
919 CHECK((nsteps % ir->nstcomm) && (ir->etc == etcANDERSENMASSIVE));
922 if (ir->etc == etcBERENDSEN)
924 sprintf(warn_buf, "The %s thermostat does not generate the correct kinetic energy distribution. You might want to consider using the %s thermostat.",
925 ETCOUPLTYPE(ir->etc), ETCOUPLTYPE(etcVRESCALE));
926 warning_note(wi, warn_buf);
929 if ((ir->etc == etcNOSEHOOVER || ETC_ANDERSEN(ir->etc))
930 && ir->epc == epcBERENDSEN)
932 sprintf(warn_buf, "Using Berendsen pressure coupling invalidates the "
933 "true ensemble for the thermostat");
934 warning(wi, warn_buf);
937 /* PRESSURE COUPLING */
938 if (ir->epc == epcISOTROPIC)
940 ir->epc = epcBERENDSEN;
941 warning_note(wi, "Old option for pressure coupling given: "
942 "changing \"Isotropic\" to \"Berendsen\"\n");
945 if (ir->epc != epcNO)
947 dt_pcoupl = ir->nstpcouple*ir->delta_t;
949 sprintf(err_buf, "tau-p must be > 0 instead of %g\n", ir->tau_p);
950 CHECK(ir->tau_p <= 0);
952 if (ir->tau_p/dt_pcoupl < pcouple_min_integration_steps(ir->epc))
954 sprintf(warn_buf, "For proper integration of the %s barostat, tau-p (%g) should be at least %d times larger than nstpcouple*dt (%g)",
955 EPCOUPLTYPE(ir->epc), ir->tau_p, pcouple_min_integration_steps(ir->epc), dt_pcoupl);
956 warning(wi, warn_buf);
959 sprintf(err_buf, "compressibility must be > 0 when using pressure"
960 " coupling %s\n", EPCOUPLTYPE(ir->epc));
961 CHECK(ir->compress[XX][XX] < 0 || ir->compress[YY][YY] < 0 ||
962 ir->compress[ZZ][ZZ] < 0 ||
963 (trace(ir->compress) == 0 && ir->compress[YY][XX] <= 0 &&
964 ir->compress[ZZ][XX] <= 0 && ir->compress[ZZ][YY] <= 0));
966 if (epcPARRINELLORAHMAN == ir->epc && opts->bGenVel)
969 "You are generating velocities so I am assuming you "
970 "are equilibrating a system. You are using "
971 "%s pressure coupling, but this can be "
972 "unstable for equilibration. If your system crashes, try "
973 "equilibrating first with Berendsen pressure coupling. If "
974 "you are not equilibrating the system, you can probably "
975 "ignore this warning.",
976 epcoupl_names[ir->epc]);
977 warning(wi, warn_buf);
985 if ((ir->epc != epcBERENDSEN) && (ir->epc != epcMTTK))
987 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.");
993 /* More checks are in triple check (grompp.c) */
995 if (ir->coulombtype == eelSWITCH)
997 sprintf(warn_buf, "coulombtype = %s is only for testing purposes and can lead to serious "
998 "artifacts, advice: use coulombtype = %s",
999 eel_names[ir->coulombtype],
1000 eel_names[eelRF_ZERO]);
1001 warning(wi, warn_buf);
1004 if (ir->epsilon_r != 1 && ir->implicit_solvent == eisGBSA)
1006 sprintf(warn_buf, "epsilon-r = %g with GB implicit solvent, will use this value for inner dielectric", ir->epsilon_r);
1007 warning_note(wi, warn_buf);
1010 if (EEL_RF(ir->coulombtype) && ir->epsilon_rf == 1 && ir->epsilon_r != 1)
1012 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);
1013 warning(wi, warn_buf);
1014 ir->epsilon_rf = ir->epsilon_r;
1015 ir->epsilon_r = 1.0;
1018 if (getenv("GALACTIC_DYNAMICS") == NULL)
1020 sprintf(err_buf, "epsilon-r must be >= 0 instead of %g\n", ir->epsilon_r);
1021 CHECK(ir->epsilon_r < 0);
1024 if (EEL_RF(ir->coulombtype))
1026 /* reaction field (at the cut-off) */
1028 if (ir->coulombtype == eelRF_ZERO)
1030 sprintf(warn_buf, "With coulombtype = %s, epsilon-rf must be 0, assuming you meant epsilon_rf=0",
1031 eel_names[ir->coulombtype]);
1032 CHECK(ir->epsilon_rf != 0);
1033 ir->epsilon_rf = 0.0;
1036 sprintf(err_buf, "epsilon-rf must be >= epsilon-r");
1037 CHECK((ir->epsilon_rf < ir->epsilon_r && ir->epsilon_rf != 0) ||
1038 (ir->epsilon_r == 0));
1039 if (ir->epsilon_rf == ir->epsilon_r)
1041 sprintf(warn_buf, "Using epsilon-rf = epsilon-r with %s does not make sense",
1042 eel_names[ir->coulombtype]);
1043 warning(wi, warn_buf);
1046 /* Allow rlist>rcoulomb for tabulated long range stuff. This just
1047 * means the interaction is zero outside rcoulomb, but it helps to
1048 * provide accurate energy conservation.
1050 if (EEL_MIGHT_BE_ZERO_AT_CUTOFF(ir->coulombtype))
1052 if (EEL_SWITCHED(ir->coulombtype))
1055 "With coulombtype = %s rcoulomb_switch must be < rcoulomb. Or, better: Use the potential modifier options!",
1056 eel_names[ir->coulombtype]);
1057 CHECK(ir->rcoulomb_switch >= ir->rcoulomb);
1060 else if (ir->coulombtype == eelCUT || EEL_RF(ir->coulombtype))
1062 if (ir->cutoff_scheme == ecutsGROUP && ir->coulomb_modifier == eintmodNONE)
1064 sprintf(err_buf, "With coulombtype = %s, rcoulomb should be >= rlist unless you use a potential modifier",
1065 eel_names[ir->coulombtype]);
1066 CHECK(ir->rlist > ir->rcoulomb);
1070 if (ir->coulombtype == eelSWITCH || ir->coulombtype == eelSHIFT ||
1071 ir->vdwtype == evdwSWITCH || ir->vdwtype == evdwSHIFT)
1074 "The switch/shift interaction settings are just for compatibility; you will get better "
1075 "performance from applying potential modifiers to your interactions!\n");
1076 warning_note(wi, warn_buf);
1079 if (EEL_FULL(ir->coulombtype))
1081 if (ir->coulombtype == eelPMESWITCH || ir->coulombtype == eelPMEUSER ||
1082 ir->coulombtype == eelPMEUSERSWITCH)
1084 sprintf(err_buf, "With coulombtype = %s, rcoulomb must be <= rlist",
1085 eel_names[ir->coulombtype]);
1086 CHECK(ir->rcoulomb > ir->rlist);
1088 else if (ir->cutoff_scheme == ecutsGROUP && ir->coulomb_modifier == eintmodNONE)
1090 if (ir->coulombtype == eelPME || ir->coulombtype == eelP3M_AD)
1093 "With coulombtype = %s (without modifier), rcoulomb must be equal to rlist,\n"
1094 "or rlistlong if nstcalclr=1. For optimal energy conservation,consider using\n"
1095 "a potential modifier.", eel_names[ir->coulombtype]);
1096 if (ir->nstcalclr == 1)
1098 CHECK(ir->rcoulomb != ir->rlist && ir->rcoulomb != ir->rlistlong);
1102 CHECK(ir->rcoulomb != ir->rlist);
1108 if (EEL_PME(ir->coulombtype))
1110 if (ir->pme_order < 3)
1112 warning_error(wi, "pme-order can not be smaller than 3");
1116 if (ir->nwall == 2 && EEL_FULL(ir->coulombtype))
1118 if (ir->ewald_geometry == eewg3D)
1120 sprintf(warn_buf, "With pbc=%s you should use ewald-geometry=%s",
1121 epbc_names[ir->ePBC], eewg_names[eewg3DC]);
1122 warning(wi, warn_buf);
1124 /* This check avoids extra pbc coding for exclusion corrections */
1125 sprintf(err_buf, "wall-ewald-zfac should be >= 2");
1126 CHECK(ir->wall_ewald_zfac < 2);
1129 if (EVDW_SWITCHED(ir->vdwtype))
1131 sprintf(err_buf, "With vdwtype = %s rvdw-switch must be < rvdw. Or, better - use a potential modifier.",
1132 evdw_names[ir->vdwtype]);
1133 CHECK(ir->rvdw_switch >= ir->rvdw);
1135 else if (ir->vdwtype == evdwCUT)
1137 if (ir->cutoff_scheme == ecutsGROUP && ir->vdw_modifier == eintmodNONE)
1139 sprintf(err_buf, "With vdwtype = %s, rvdw must be >= rlist unless you use a potential modifier", evdw_names[ir->vdwtype]);
1140 CHECK(ir->rlist > ir->rvdw);
1143 if (ir->cutoff_scheme == ecutsGROUP)
1145 if (EEL_IS_ZERO_AT_CUTOFF(ir->coulombtype)
1146 && (ir->rlistlong <= ir->rcoulomb))
1148 sprintf(warn_buf, "For energy conservation with switch/shift potentials, %s should be 0.1 to 0.3 nm larger than rcoulomb.",
1149 IR_TWINRANGE(*ir) ? "rlistlong" : "rlist");
1150 warning_note(wi, warn_buf);
1152 if (EVDW_SWITCHED(ir->vdwtype) && (ir->rlistlong <= ir->rvdw))
1154 sprintf(warn_buf, "For energy conservation with switch/shift potentials, %s should be 0.1 to 0.3 nm larger than rvdw.",
1155 IR_TWINRANGE(*ir) ? "rlistlong" : "rlist");
1156 warning_note(wi, warn_buf);
1160 if (ir->vdwtype == evdwUSER && ir->eDispCorr != edispcNO)
1162 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.");
1165 if (ir->nstlist == -1)
1167 sprintf(err_buf, "With nstlist=-1 rvdw and rcoulomb should be smaller than rlist to account for diffusion and possibly charge-group radii");
1168 CHECK(ir->rvdw >= ir->rlist || ir->rcoulomb >= ir->rlist);
1170 sprintf(err_buf, "nstlist can not be smaller than -1");
1171 CHECK(ir->nstlist < -1);
1173 if (ir->eI == eiLBFGS && (ir->coulombtype == eelCUT || ir->vdwtype == evdwCUT)
1176 warning(wi, "For efficient BFGS minimization, use switch/shift/pme instead of cut-off.");
1179 if (ir->eI == eiLBFGS && ir->nbfgscorr <= 0)
1181 warning(wi, "Using L-BFGS with nbfgscorr<=0 just gets you steepest descent.");
1184 /* ENERGY CONSERVATION */
1185 if (ir_NVE(ir) && ir->cutoff_scheme == ecutsGROUP)
1187 if (!EVDW_MIGHT_BE_ZERO_AT_CUTOFF(ir->vdwtype) && ir->rvdw > 0 && ir->vdw_modifier == eintmodNONE)
1189 sprintf(warn_buf, "You are using a cut-off for VdW interactions with NVE, for good energy conservation use vdwtype = %s (possibly with DispCorr)",
1190 evdw_names[evdwSHIFT]);
1191 warning_note(wi, warn_buf);
1193 if (!EEL_MIGHT_BE_ZERO_AT_CUTOFF(ir->coulombtype) && ir->rcoulomb > 0 && ir->coulomb_modifier == eintmodNONE)
1195 sprintf(warn_buf, "You are using a cut-off for electrostatics with NVE, for good energy conservation use coulombtype = %s or %s",
1196 eel_names[eelPMESWITCH], eel_names[eelRF_ZERO]);
1197 warning_note(wi, warn_buf);
1201 /* IMPLICIT SOLVENT */
1202 if (ir->coulombtype == eelGB_NOTUSED)
1204 ir->coulombtype = eelCUT;
1205 ir->implicit_solvent = eisGBSA;
1206 fprintf(stderr, "Note: Old option for generalized born electrostatics given:\n"
1207 "Changing coulombtype from \"generalized-born\" to \"cut-off\" and instead\n"
1208 "setting implicit-solvent value to \"GBSA\" in input section.\n");
1211 if (ir->sa_algorithm == esaSTILL)
1213 sprintf(err_buf, "Still SA algorithm not available yet, use %s or %s instead\n", esa_names[esaAPPROX], esa_names[esaNO]);
1214 CHECK(ir->sa_algorithm == esaSTILL);
1217 if (ir->implicit_solvent == eisGBSA)
1219 sprintf(err_buf, "With GBSA implicit solvent, rgbradii must be equal to rlist.");
1220 CHECK(ir->rgbradii != ir->rlist);
1222 if (ir->coulombtype != eelCUT)
1224 sprintf(err_buf, "With GBSA, coulombtype must be equal to %s\n", eel_names[eelCUT]);
1225 CHECK(ir->coulombtype != eelCUT);
1227 if (ir->vdwtype != evdwCUT)
1229 sprintf(err_buf, "With GBSA, vdw-type must be equal to %s\n", evdw_names[evdwCUT]);
1230 CHECK(ir->vdwtype != evdwCUT);
1232 if (ir->nstgbradii < 1)
1234 sprintf(warn_buf, "Using GBSA with nstgbradii<1, setting nstgbradii=1");
1235 warning_note(wi, warn_buf);
1238 if (ir->sa_algorithm == esaNO)
1240 sprintf(warn_buf, "No SA (non-polar) calculation requested together with GB. Are you sure this is what you want?\n");
1241 warning_note(wi, warn_buf);
1243 if (ir->sa_surface_tension < 0 && ir->sa_algorithm != esaNO)
1245 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");
1246 warning_note(wi, warn_buf);
1248 if (ir->gb_algorithm == egbSTILL)
1250 ir->sa_surface_tension = 0.0049 * CAL2JOULE * 100;
1254 ir->sa_surface_tension = 0.0054 * CAL2JOULE * 100;
1257 if (ir->sa_surface_tension == 0 && ir->sa_algorithm != esaNO)
1259 sprintf(err_buf, "Surface tension set to 0 while SA-calculation requested\n");
1260 CHECK(ir->sa_surface_tension == 0 && ir->sa_algorithm != esaNO);
1267 if (ir->cutoff_scheme != ecutsGROUP)
1269 warning_error(wi, "AdresS simulation supports only cutoff-scheme=group");
1273 warning_error(wi, "AdresS simulation supports only stochastic dynamics");
1275 if (ir->epc != epcNO)
1277 warning_error(wi, "AdresS simulation does not support pressure coupling");
1279 if (EEL_FULL(ir->coulombtype))
1281 warning_error(wi, "AdresS simulation does not support long-range electrostatics");
1286 /* count the number of text elemets separated by whitespace in a string.
1287 str = the input string
1288 maxptr = the maximum number of allowed elements
1289 ptr = the output array of pointers to the first character of each element
1290 returns: the number of elements. */
1291 int str_nelem(const char *str, int maxptr, char *ptr[])
1296 copy0 = strdup(str);
1299 while (*copy != '\0')
1303 gmx_fatal(FARGS, "Too many groups on line: '%s' (max is %d)",
1311 while ((*copy != '\0') && !isspace(*copy))
1330 /* interpret a number of doubles from a string and put them in an array,
1331 after allocating space for them.
1332 str = the input string
1333 n = the (pre-allocated) number of doubles read
1334 r = the output array of doubles. */
1335 static void parse_n_real(char *str, int *n, real **r)
1340 *n = str_nelem(str, MAXPTR, ptr);
1343 for (i = 0; i < *n; i++)
1345 (*r)[i] = strtod(ptr[i], NULL);
1349 static void do_fep_params(t_inputrec *ir, char fep_lambda[][STRLEN], char weights[STRLEN])
1352 int i, j, max_n_lambda, nweights, nfep[efptNR];
1353 t_lambda *fep = ir->fepvals;
1354 t_expanded *expand = ir->expandedvals;
1355 real **count_fep_lambdas;
1356 gmx_bool bOneLambda = TRUE;
1358 snew(count_fep_lambdas, efptNR);
1360 /* FEP input processing */
1361 /* first, identify the number of lambda values for each type.
1362 All that are nonzero must have the same number */
1364 for (i = 0; i < efptNR; i++)
1366 parse_n_real(fep_lambda[i], &(nfep[i]), &(count_fep_lambdas[i]));
1369 /* now, determine the number of components. All must be either zero, or equal. */
1372 for (i = 0; i < efptNR; i++)
1374 if (nfep[i] > max_n_lambda)
1376 max_n_lambda = nfep[i]; /* here's a nonzero one. All of them
1377 must have the same number if its not zero.*/
1382 for (i = 0; i < efptNR; i++)
1386 ir->fepvals->separate_dvdl[i] = FALSE;
1388 else if (nfep[i] == max_n_lambda)
1390 if (i != efptTEMPERATURE) /* we treat this differently -- not really a reason to compute the derivative with
1391 respect to the temperature currently */
1393 ir->fepvals->separate_dvdl[i] = TRUE;
1398 gmx_fatal(FARGS, "Number of lambdas (%d) for FEP type %s not equal to number of other types (%d)",
1399 nfep[i], efpt_names[i], max_n_lambda);
1402 /* we don't print out dhdl if the temperature is changing, since we can't correctly define dhdl in this case */
1403 ir->fepvals->separate_dvdl[efptTEMPERATURE] = FALSE;
1405 /* the number of lambdas is the number we've read in, which is either zero
1406 or the same for all */
1407 fep->n_lambda = max_n_lambda;
1409 /* allocate space for the array of lambda values */
1410 snew(fep->all_lambda, efptNR);
1411 /* if init_lambda is defined, we need to set lambda */
1412 if ((fep->init_lambda > 0) && (fep->n_lambda == 0))
1414 ir->fepvals->separate_dvdl[efptFEP] = TRUE;
1416 /* otherwise allocate the space for all of the lambdas, and transfer the data */
1417 for (i = 0; i < efptNR; i++)
1419 snew(fep->all_lambda[i], fep->n_lambda);
1420 if (nfep[i] > 0) /* if it's zero, then the count_fep_lambda arrays
1423 for (j = 0; j < fep->n_lambda; j++)
1425 fep->all_lambda[i][j] = (double)count_fep_lambdas[i][j];
1427 sfree(count_fep_lambdas[i]);
1430 sfree(count_fep_lambdas);
1432 /* "fep-vals" is either zero or the full number. If zero, we'll need to define fep-lambdas for internal
1433 bookkeeping -- for now, init_lambda */
1435 if ((nfep[efptFEP] == 0) && (fep->init_lambda >= 0))
1437 for (i = 0; i < fep->n_lambda; i++)
1439 fep->all_lambda[efptFEP][i] = fep->init_lambda;
1443 /* check to see if only a single component lambda is defined, and soft core is defined.
1444 In this case, turn on coulomb soft core */
1446 if (max_n_lambda == 0)
1452 for (i = 0; i < efptNR; i++)
1454 if ((nfep[i] != 0) && (i != efptFEP))
1460 if ((bOneLambda) && (fep->sc_alpha > 0))
1462 fep->bScCoul = TRUE;
1465 /* Fill in the others with the efptFEP if they are not explicitly
1466 specified (i.e. nfep[i] == 0). This means if fep is not defined,
1467 they are all zero. */
1469 for (i = 0; i < efptNR; i++)
1471 if ((nfep[i] == 0) && (i != efptFEP))
1473 for (j = 0; j < fep->n_lambda; j++)
1475 fep->all_lambda[i][j] = fep->all_lambda[efptFEP][j];
1481 /* make it easier if sc_r_power = 48 by increasing it to the 4th power, to be in the right scale. */
1482 if (fep->sc_r_power == 48)
1484 if (fep->sc_alpha > 0.1)
1486 gmx_fatal(FARGS, "sc_alpha (%f) for sc_r_power = 48 should usually be between 0.001 and 0.004", fep->sc_alpha);
1490 expand = ir->expandedvals;
1491 /* now read in the weights */
1492 parse_n_real(weights, &nweights, &(expand->init_lambda_weights));
1495 expand->bInit_weights = FALSE;
1496 snew(expand->init_lambda_weights, fep->n_lambda); /* initialize to zero */
1498 else if (nweights != fep->n_lambda)
1500 gmx_fatal(FARGS, "Number of weights (%d) is not equal to number of lambda values (%d)",
1501 nweights, fep->n_lambda);
1505 expand->bInit_weights = TRUE;
1507 if ((expand->nstexpanded < 0) && (ir->efep != efepNO))
1509 expand->nstexpanded = fep->nstdhdl;
1510 /* if you don't specify nstexpanded when doing expanded ensemble free energy calcs, it is set to nstdhdl */
1512 if ((expand->nstexpanded < 0) && ir->bSimTemp)
1514 expand->nstexpanded = 2*(int)(ir->opts.tau_t[0]/ir->delta_t);
1515 /* if you don't specify nstexpanded when doing expanded ensemble simulated tempering, it is set to
1516 2*tau_t just to be careful so it's not to frequent */
1521 static void do_simtemp_params(t_inputrec *ir)
1524 snew(ir->simtempvals->temperatures, ir->fepvals->n_lambda);
1525 GetSimTemps(ir->fepvals->n_lambda, ir->simtempvals, ir->fepvals->all_lambda[efptTEMPERATURE]);
1530 static void do_wall_params(t_inputrec *ir,
1531 char *wall_atomtype, char *wall_density,
1535 char *names[MAXPTR];
1538 opts->wall_atomtype[0] = NULL;
1539 opts->wall_atomtype[1] = NULL;
1541 ir->wall_atomtype[0] = -1;
1542 ir->wall_atomtype[1] = -1;
1543 ir->wall_density[0] = 0;
1544 ir->wall_density[1] = 0;
1548 nstr = str_nelem(wall_atomtype, MAXPTR, names);
1549 if (nstr != ir->nwall)
1551 gmx_fatal(FARGS, "Expected %d elements for wall_atomtype, found %d",
1554 for (i = 0; i < ir->nwall; i++)
1556 opts->wall_atomtype[i] = strdup(names[i]);
1559 if (ir->wall_type == ewt93 || ir->wall_type == ewt104)
1561 nstr = str_nelem(wall_density, MAXPTR, names);
1562 if (nstr != ir->nwall)
1564 gmx_fatal(FARGS, "Expected %d elements for wall-density, found %d", ir->nwall, nstr);
1566 for (i = 0; i < ir->nwall; i++)
1568 sscanf(names[i], "%lf", &dbl);
1571 gmx_fatal(FARGS, "wall-density[%d] = %f\n", i, dbl);
1573 ir->wall_density[i] = dbl;
1579 static void add_wall_energrps(gmx_groups_t *groups, int nwall, t_symtab *symtab)
1587 srenew(groups->grpname, groups->ngrpname+nwall);
1588 grps = &(groups->grps[egcENER]);
1589 srenew(grps->nm_ind, grps->nr+nwall);
1590 for (i = 0; i < nwall; i++)
1592 sprintf(str, "wall%d", i);
1593 groups->grpname[groups->ngrpname] = put_symtab(symtab, str);
1594 grps->nm_ind[grps->nr++] = groups->ngrpname++;
1599 void read_expandedparams(int *ninp_p, t_inpfile **inp_p,
1600 t_expanded *expand, warninp_t wi)
1602 int ninp, nerror = 0;
1608 /* read expanded ensemble parameters */
1609 CCTYPE ("expanded ensemble variables");
1610 ITYPE ("nstexpanded", expand->nstexpanded, -1);
1611 EETYPE("lmc-stats", expand->elamstats, elamstats_names);
1612 EETYPE("lmc-move", expand->elmcmove, elmcmove_names);
1613 EETYPE("lmc-weights-equil", expand->elmceq, elmceq_names);
1614 ITYPE ("weight-equil-number-all-lambda", expand->equil_n_at_lam, -1);
1615 ITYPE ("weight-equil-number-samples", expand->equil_samples, -1);
1616 ITYPE ("weight-equil-number-steps", expand->equil_steps, -1);
1617 RTYPE ("weight-equil-wl-delta", expand->equil_wl_delta, -1);
1618 RTYPE ("weight-equil-count-ratio", expand->equil_ratio, -1);
1619 CCTYPE("Seed for Monte Carlo in lambda space");
1620 ITYPE ("lmc-seed", expand->lmc_seed, -1);
1621 RTYPE ("mc-temperature", expand->mc_temp, -1);
1622 ITYPE ("lmc-repeats", expand->lmc_repeats, 1);
1623 ITYPE ("lmc-gibbsdelta", expand->gibbsdeltalam, -1);
1624 ITYPE ("lmc-forced-nstart", expand->lmc_forced_nstart, 0);
1625 EETYPE("symmetrized-transition-matrix", expand->bSymmetrizedTMatrix, yesno_names);
1626 ITYPE("nst-transition-matrix", expand->nstTij, -1);
1627 ITYPE ("mininum-var-min", expand->minvarmin, 100); /*default is reasonable */
1628 ITYPE ("weight-c-range", expand->c_range, 0); /* default is just C=0 */
1629 RTYPE ("wl-scale", expand->wl_scale, 0.8);
1630 RTYPE ("wl-ratio", expand->wl_ratio, 0.8);
1631 RTYPE ("init-wl-delta", expand->init_wl_delta, 1.0);
1632 EETYPE("wl-oneovert", expand->bWLoneovert, yesno_names);
1640 void get_ir(const char *mdparin, const char *mdparout,
1641 t_inputrec *ir, t_gromppopts *opts,
1645 double dumdub[2][6];
1649 char warn_buf[STRLEN];
1650 t_lambda *fep = ir->fepvals;
1651 t_expanded *expand = ir->expandedvals;
1653 inp = read_inpfile(mdparin, &ninp, NULL, wi);
1655 snew(dumstr[0], STRLEN);
1656 snew(dumstr[1], STRLEN);
1658 /* remove the following deprecated commands */
1661 REM_TYPE("domain-decomposition");
1662 REM_TYPE("andersen-seed");
1664 REM_TYPE("dihre-fc");
1665 REM_TYPE("dihre-tau");
1666 REM_TYPE("nstdihreout");
1667 REM_TYPE("nstcheckpoint");
1669 /* replace the following commands with the clearer new versions*/
1670 REPL_TYPE("unconstrained-start", "continuation");
1671 REPL_TYPE("foreign-lambda", "fep-lambdas");
1673 CCTYPE ("VARIOUS PREPROCESSING OPTIONS");
1674 CTYPE ("Preprocessor information: use cpp syntax.");
1675 CTYPE ("e.g.: -I/home/joe/doe -I/home/mary/roe");
1676 STYPE ("include", opts->include, NULL);
1677 CTYPE ("e.g.: -DPOSRES -DFLEXIBLE (note these variable names are case sensitive)");
1678 STYPE ("define", opts->define, NULL);
1680 CCTYPE ("RUN CONTROL PARAMETERS");
1681 EETYPE("integrator", ir->eI, ei_names);
1682 CTYPE ("Start time and timestep in ps");
1683 RTYPE ("tinit", ir->init_t, 0.0);
1684 RTYPE ("dt", ir->delta_t, 0.001);
1685 STEPTYPE ("nsteps", ir->nsteps, 0);
1686 CTYPE ("For exact run continuation or redoing part of a run");
1687 STEPTYPE ("init-step", ir->init_step, 0);
1688 CTYPE ("Part index is updated automatically on checkpointing (keeps files separate)");
1689 ITYPE ("simulation-part", ir->simulation_part, 1);
1690 CTYPE ("mode for center of mass motion removal");
1691 EETYPE("comm-mode", ir->comm_mode, ecm_names);
1692 CTYPE ("number of steps for center of mass motion removal");
1693 ITYPE ("nstcomm", ir->nstcomm, 100);
1694 CTYPE ("group(s) for center of mass motion removal");
1695 STYPE ("comm-grps", vcm, NULL);
1697 CCTYPE ("LANGEVIN DYNAMICS OPTIONS");
1698 CTYPE ("Friction coefficient (amu/ps) and random seed");
1699 RTYPE ("bd-fric", ir->bd_fric, 0.0);
1700 ITYPE ("ld-seed", ir->ld_seed, 1993);
1703 CCTYPE ("ENERGY MINIMIZATION OPTIONS");
1704 CTYPE ("Force tolerance and initial step-size");
1705 RTYPE ("emtol", ir->em_tol, 10.0);
1706 RTYPE ("emstep", ir->em_stepsize, 0.01);
1707 CTYPE ("Max number of iterations in relax-shells");
1708 ITYPE ("niter", ir->niter, 20);
1709 CTYPE ("Step size (ps^2) for minimization of flexible constraints");
1710 RTYPE ("fcstep", ir->fc_stepsize, 0);
1711 CTYPE ("Frequency of steepest descents steps when doing CG");
1712 ITYPE ("nstcgsteep", ir->nstcgsteep, 1000);
1713 ITYPE ("nbfgscorr", ir->nbfgscorr, 10);
1715 CCTYPE ("TEST PARTICLE INSERTION OPTIONS");
1716 RTYPE ("rtpi", ir->rtpi, 0.05);
1718 /* Output options */
1719 CCTYPE ("OUTPUT CONTROL OPTIONS");
1720 CTYPE ("Output frequency for coords (x), velocities (v) and forces (f)");
1721 ITYPE ("nstxout", ir->nstxout, 0);
1722 ITYPE ("nstvout", ir->nstvout, 0);
1723 ITYPE ("nstfout", ir->nstfout, 0);
1724 ir->nstcheckpoint = 1000;
1725 CTYPE ("Output frequency for energies to log file and energy file");
1726 ITYPE ("nstlog", ir->nstlog, 1000);
1727 ITYPE ("nstcalcenergy", ir->nstcalcenergy, 100);
1728 ITYPE ("nstenergy", ir->nstenergy, 1000);
1729 CTYPE ("Output frequency and precision for .xtc file");
1730 ITYPE ("nstxtcout", ir->nstxtcout, 0);
1731 RTYPE ("xtc-precision", ir->xtcprec, 1000.0);
1732 CTYPE ("This selects the subset of atoms for the .xtc file. You can");
1733 CTYPE ("select multiple groups. By default all atoms will be written.");
1734 STYPE ("xtc-grps", xtc_grps, NULL);
1735 CTYPE ("Selection of energy groups");
1736 STYPE ("energygrps", energy, NULL);
1738 /* Neighbor searching */
1739 CCTYPE ("NEIGHBORSEARCHING PARAMETERS");
1740 CTYPE ("cut-off scheme (group: using charge groups, Verlet: particle based cut-offs)");
1741 EETYPE("cutoff-scheme", ir->cutoff_scheme, ecutscheme_names);
1742 CTYPE ("nblist update frequency");
1743 ITYPE ("nstlist", ir->nstlist, 10);
1744 CTYPE ("ns algorithm (simple or grid)");
1745 EETYPE("ns-type", ir->ns_type, ens_names);
1746 /* set ndelta to the optimal value of 2 */
1748 CTYPE ("Periodic boundary conditions: xyz, no, xy");
1749 EETYPE("pbc", ir->ePBC, epbc_names);
1750 EETYPE("periodic-molecules", ir->bPeriodicMols, yesno_names);
1751 CTYPE ("Allowed energy drift due to the Verlet buffer in kJ/mol/ps per atom,");
1752 CTYPE ("a value of -1 means: use rlist");
1753 RTYPE("verlet-buffer-drift", ir->verletbuf_drift, 0.005);
1754 CTYPE ("nblist cut-off");
1755 RTYPE ("rlist", ir->rlist, 1.0);
1756 CTYPE ("long-range cut-off for switched potentials");
1757 RTYPE ("rlistlong", ir->rlistlong, -1);
1758 ITYPE ("nstcalclr", ir->nstcalclr, -1);
1760 /* Electrostatics */
1761 CCTYPE ("OPTIONS FOR ELECTROSTATICS AND VDW");
1762 CTYPE ("Method for doing electrostatics");
1763 EETYPE("coulombtype", ir->coulombtype, eel_names);
1764 EETYPE("coulomb-modifier", ir->coulomb_modifier, eintmod_names);
1765 CTYPE ("cut-off lengths");
1766 RTYPE ("rcoulomb-switch", ir->rcoulomb_switch, 0.0);
1767 RTYPE ("rcoulomb", ir->rcoulomb, 1.0);
1768 CTYPE ("Relative dielectric constant for the medium and the reaction field");
1769 RTYPE ("epsilon-r", ir->epsilon_r, 1.0);
1770 RTYPE ("epsilon-rf", ir->epsilon_rf, 0.0);
1771 CTYPE ("Method for doing Van der Waals");
1772 EETYPE("vdw-type", ir->vdwtype, evdw_names);
1773 EETYPE("vdw-modifier", ir->vdw_modifier, eintmod_names);
1774 CTYPE ("cut-off lengths");
1775 RTYPE ("rvdw-switch", ir->rvdw_switch, 0.0);
1776 RTYPE ("rvdw", ir->rvdw, 1.0);
1777 CTYPE ("Apply long range dispersion corrections for Energy and Pressure");
1778 EETYPE("DispCorr", ir->eDispCorr, edispc_names);
1779 CTYPE ("Extension of the potential lookup tables beyond the cut-off");
1780 RTYPE ("table-extension", ir->tabext, 1.0);
1781 CTYPE ("Separate tables between energy group pairs");
1782 STYPE ("energygrp-table", egptable, NULL);
1783 CTYPE ("Spacing for the PME/PPPM FFT grid");
1784 RTYPE ("fourierspacing", ir->fourier_spacing, 0.12);
1785 CTYPE ("FFT grid size, when a value is 0 fourierspacing will be used");
1786 ITYPE ("fourier-nx", ir->nkx, 0);
1787 ITYPE ("fourier-ny", ir->nky, 0);
1788 ITYPE ("fourier-nz", ir->nkz, 0);
1789 CTYPE ("EWALD/PME/PPPM parameters");
1790 ITYPE ("pme-order", ir->pme_order, 4);
1791 RTYPE ("ewald-rtol", ir->ewald_rtol, 0.00001);
1792 EETYPE("ewald-geometry", ir->ewald_geometry, eewg_names);
1793 RTYPE ("epsilon-surface", ir->epsilon_surface, 0.0);
1794 EETYPE("optimize-fft", ir->bOptFFT, yesno_names);
1796 CCTYPE("IMPLICIT SOLVENT ALGORITHM");
1797 EETYPE("implicit-solvent", ir->implicit_solvent, eis_names);
1799 CCTYPE ("GENERALIZED BORN ELECTROSTATICS");
1800 CTYPE ("Algorithm for calculating Born radii");
1801 EETYPE("gb-algorithm", ir->gb_algorithm, egb_names);
1802 CTYPE ("Frequency of calculating the Born radii inside rlist");
1803 ITYPE ("nstgbradii", ir->nstgbradii, 1);
1804 CTYPE ("Cutoff for Born radii calculation; the contribution from atoms");
1805 CTYPE ("between rlist and rgbradii is updated every nstlist steps");
1806 RTYPE ("rgbradii", ir->rgbradii, 1.0);
1807 CTYPE ("Dielectric coefficient of the implicit solvent");
1808 RTYPE ("gb-epsilon-solvent", ir->gb_epsilon_solvent, 80.0);
1809 CTYPE ("Salt concentration in M for Generalized Born models");
1810 RTYPE ("gb-saltconc", ir->gb_saltconc, 0.0);
1811 CTYPE ("Scaling factors used in the OBC GB model. Default values are OBC(II)");
1812 RTYPE ("gb-obc-alpha", ir->gb_obc_alpha, 1.0);
1813 RTYPE ("gb-obc-beta", ir->gb_obc_beta, 0.8);
1814 RTYPE ("gb-obc-gamma", ir->gb_obc_gamma, 4.85);
1815 RTYPE ("gb-dielectric-offset", ir->gb_dielectric_offset, 0.009);
1816 EETYPE("sa-algorithm", ir->sa_algorithm, esa_names);
1817 CTYPE ("Surface tension (kJ/mol/nm^2) for the SA (nonpolar surface) part of GBSA");
1818 CTYPE ("The value -1 will set default value for Still/HCT/OBC GB-models.");
1819 RTYPE ("sa-surface-tension", ir->sa_surface_tension, -1);
1821 /* Coupling stuff */
1822 CCTYPE ("OPTIONS FOR WEAK COUPLING ALGORITHMS");
1823 CTYPE ("Temperature coupling");
1824 EETYPE("tcoupl", ir->etc, etcoupl_names);
1825 ITYPE ("nsttcouple", ir->nsttcouple, -1);
1826 ITYPE("nh-chain-length", ir->opts.nhchainlength, NHCHAINLENGTH);
1827 EETYPE("print-nose-hoover-chain-variables", ir->bPrintNHChains, yesno_names);
1828 CTYPE ("Groups to couple separately");
1829 STYPE ("tc-grps", tcgrps, NULL);
1830 CTYPE ("Time constant (ps) and reference temperature (K)");
1831 STYPE ("tau-t", tau_t, NULL);
1832 STYPE ("ref-t", ref_t, NULL);
1833 CTYPE ("pressure coupling");
1834 EETYPE("pcoupl", ir->epc, epcoupl_names);
1835 EETYPE("pcoupltype", ir->epct, epcoupltype_names);
1836 ITYPE ("nstpcouple", ir->nstpcouple, -1);
1837 CTYPE ("Time constant (ps), compressibility (1/bar) and reference P (bar)");
1838 RTYPE ("tau-p", ir->tau_p, 1.0);
1839 STYPE ("compressibility", dumstr[0], NULL);
1840 STYPE ("ref-p", dumstr[1], NULL);
1841 CTYPE ("Scaling of reference coordinates, No, All or COM");
1842 EETYPE ("refcoord-scaling", ir->refcoord_scaling, erefscaling_names);
1845 CCTYPE ("OPTIONS FOR QMMM calculations");
1846 EETYPE("QMMM", ir->bQMMM, yesno_names);
1847 CTYPE ("Groups treated Quantum Mechanically");
1848 STYPE ("QMMM-grps", QMMM, NULL);
1849 CTYPE ("QM method");
1850 STYPE("QMmethod", QMmethod, NULL);
1851 CTYPE ("QMMM scheme");
1852 EETYPE("QMMMscheme", ir->QMMMscheme, eQMMMscheme_names);
1853 CTYPE ("QM basisset");
1854 STYPE("QMbasis", QMbasis, NULL);
1855 CTYPE ("QM charge");
1856 STYPE ("QMcharge", QMcharge, NULL);
1857 CTYPE ("QM multiplicity");
1858 STYPE ("QMmult", QMmult, NULL);
1859 CTYPE ("Surface Hopping");
1860 STYPE ("SH", bSH, NULL);
1861 CTYPE ("CAS space options");
1862 STYPE ("CASorbitals", CASorbitals, NULL);
1863 STYPE ("CASelectrons", CASelectrons, NULL);
1864 STYPE ("SAon", SAon, NULL);
1865 STYPE ("SAoff", SAoff, NULL);
1866 STYPE ("SAsteps", SAsteps, NULL);
1867 CTYPE ("Scale factor for MM charges");
1868 RTYPE ("MMChargeScaleFactor", ir->scalefactor, 1.0);
1869 CTYPE ("Optimization of QM subsystem");
1870 STYPE ("bOPT", bOPT, NULL);
1871 STYPE ("bTS", bTS, NULL);
1873 /* Simulated annealing */
1874 CCTYPE("SIMULATED ANNEALING");
1875 CTYPE ("Type of annealing for each temperature group (no/single/periodic)");
1876 STYPE ("annealing", anneal, NULL);
1877 CTYPE ("Number of time points to use for specifying annealing in each group");
1878 STYPE ("annealing-npoints", anneal_npoints, NULL);
1879 CTYPE ("List of times at the annealing points for each group");
1880 STYPE ("annealing-time", anneal_time, NULL);
1881 CTYPE ("Temp. at each annealing point, for each group.");
1882 STYPE ("annealing-temp", anneal_temp, NULL);
1885 CCTYPE ("GENERATE VELOCITIES FOR STARTUP RUN");
1886 EETYPE("gen-vel", opts->bGenVel, yesno_names);
1887 RTYPE ("gen-temp", opts->tempi, 300.0);
1888 ITYPE ("gen-seed", opts->seed, 173529);
1891 CCTYPE ("OPTIONS FOR BONDS");
1892 EETYPE("constraints", opts->nshake, constraints);
1893 CTYPE ("Type of constraint algorithm");
1894 EETYPE("constraint-algorithm", ir->eConstrAlg, econstr_names);
1895 CTYPE ("Do not constrain the start configuration");
1896 EETYPE("continuation", ir->bContinuation, yesno_names);
1897 CTYPE ("Use successive overrelaxation to reduce the number of shake iterations");
1898 EETYPE("Shake-SOR", ir->bShakeSOR, yesno_names);
1899 CTYPE ("Relative tolerance of shake");
1900 RTYPE ("shake-tol", ir->shake_tol, 0.0001);
1901 CTYPE ("Highest order in the expansion of the constraint coupling matrix");
1902 ITYPE ("lincs-order", ir->nProjOrder, 4);
1903 CTYPE ("Number of iterations in the final step of LINCS. 1 is fine for");
1904 CTYPE ("normal simulations, but use 2 to conserve energy in NVE runs.");
1905 CTYPE ("For energy minimization with constraints it should be 4 to 8.");
1906 ITYPE ("lincs-iter", ir->nLincsIter, 1);
1907 CTYPE ("Lincs will write a warning to the stderr if in one step a bond");
1908 CTYPE ("rotates over more degrees than");
1909 RTYPE ("lincs-warnangle", ir->LincsWarnAngle, 30.0);
1910 CTYPE ("Convert harmonic bonds to morse potentials");
1911 EETYPE("morse", opts->bMorse, yesno_names);
1913 /* Energy group exclusions */
1914 CCTYPE ("ENERGY GROUP EXCLUSIONS");
1915 CTYPE ("Pairs of energy groups for which all non-bonded interactions are excluded");
1916 STYPE ("energygrp-excl", egpexcl, NULL);
1920 CTYPE ("Number of walls, type, atom types, densities and box-z scale factor for Ewald");
1921 ITYPE ("nwall", ir->nwall, 0);
1922 EETYPE("wall-type", ir->wall_type, ewt_names);
1923 RTYPE ("wall-r-linpot", ir->wall_r_linpot, -1);
1924 STYPE ("wall-atomtype", wall_atomtype, NULL);
1925 STYPE ("wall-density", wall_density, NULL);
1926 RTYPE ("wall-ewald-zfac", ir->wall_ewald_zfac, 3);
1929 CCTYPE("COM PULLING");
1930 CTYPE("Pull type: no, umbrella, constraint or constant-force");
1931 EETYPE("pull", ir->ePull, epull_names);
1932 if (ir->ePull != epullNO)
1935 pull_grp = read_pullparams(&ninp, &inp, ir->pull, &opts->pull_start, wi);
1938 /* Enforced rotation */
1939 CCTYPE("ENFORCED ROTATION");
1940 CTYPE("Enforced rotation: No or Yes");
1941 EETYPE("rotation", ir->bRot, yesno_names);
1945 rot_grp = read_rotparams(&ninp, &inp, ir->rot, wi);
1949 CCTYPE("NMR refinement stuff");
1950 CTYPE ("Distance restraints type: No, Simple or Ensemble");
1951 EETYPE("disre", ir->eDisre, edisre_names);
1952 CTYPE ("Force weighting of pairs in one distance restraint: Conservative or Equal");
1953 EETYPE("disre-weighting", ir->eDisreWeighting, edisreweighting_names);
1954 CTYPE ("Use sqrt of the time averaged times the instantaneous violation");
1955 EETYPE("disre-mixed", ir->bDisreMixed, yesno_names);
1956 RTYPE ("disre-fc", ir->dr_fc, 1000.0);
1957 RTYPE ("disre-tau", ir->dr_tau, 0.0);
1958 CTYPE ("Output frequency for pair distances to energy file");
1959 ITYPE ("nstdisreout", ir->nstdisreout, 100);
1960 CTYPE ("Orientation restraints: No or Yes");
1961 EETYPE("orire", opts->bOrire, yesno_names);
1962 CTYPE ("Orientation restraints force constant and tau for time averaging");
1963 RTYPE ("orire-fc", ir->orires_fc, 0.0);
1964 RTYPE ("orire-tau", ir->orires_tau, 0.0);
1965 STYPE ("orire-fitgrp", orirefitgrp, NULL);
1966 CTYPE ("Output frequency for trace(SD) and S to energy file");
1967 ITYPE ("nstorireout", ir->nstorireout, 100);
1969 /* free energy variables */
1970 CCTYPE ("Free energy variables");
1971 EETYPE("free-energy", ir->efep, efep_names);
1972 STYPE ("couple-moltype", couple_moltype, NULL);
1973 EETYPE("couple-lambda0", opts->couple_lam0, couple_lam);
1974 EETYPE("couple-lambda1", opts->couple_lam1, couple_lam);
1975 EETYPE("couple-intramol", opts->bCoupleIntra, yesno_names);
1977 RTYPE ("init-lambda", fep->init_lambda, -1); /* start with -1 so
1979 it was not entered */
1980 ITYPE ("init-lambda-state", fep->init_fep_state, -1);
1981 RTYPE ("delta-lambda", fep->delta_lambda, 0.0);
1982 ITYPE ("nstdhdl", fep->nstdhdl, 50);
1983 STYPE ("fep-lambdas", fep_lambda[efptFEP], NULL);
1984 STYPE ("mass-lambdas", fep_lambda[efptMASS], NULL);
1985 STYPE ("coul-lambdas", fep_lambda[efptCOUL], NULL);
1986 STYPE ("vdw-lambdas", fep_lambda[efptVDW], NULL);
1987 STYPE ("bonded-lambdas", fep_lambda[efptBONDED], NULL);
1988 STYPE ("restraint-lambdas", fep_lambda[efptRESTRAINT], NULL);
1989 STYPE ("temperature-lambdas", fep_lambda[efptTEMPERATURE], NULL);
1990 ITYPE ("calc-lambda-neighbors", fep->lambda_neighbors, 1);
1991 STYPE ("init-lambda-weights", lambda_weights, NULL);
1992 EETYPE("dhdl-print-energy", fep->bPrintEnergy, yesno_names);
1993 RTYPE ("sc-alpha", fep->sc_alpha, 0.0);
1994 ITYPE ("sc-power", fep->sc_power, 1);
1995 RTYPE ("sc-r-power", fep->sc_r_power, 6.0);
1996 RTYPE ("sc-sigma", fep->sc_sigma, 0.3);
1997 EETYPE("sc-coul", fep->bScCoul, yesno_names);
1998 ITYPE ("dh_hist_size", fep->dh_hist_size, 0);
1999 RTYPE ("dh_hist_spacing", fep->dh_hist_spacing, 0.1);
2000 EETYPE("separate-dhdl-file", fep->separate_dhdl_file,
2001 separate_dhdl_file_names);
2002 EETYPE("dhdl-derivatives", fep->dhdl_derivatives, dhdl_derivatives_names);
2003 ITYPE ("dh_hist_size", fep->dh_hist_size, 0);
2004 RTYPE ("dh_hist_spacing", fep->dh_hist_spacing, 0.1);
2006 /* Non-equilibrium MD stuff */
2007 CCTYPE("Non-equilibrium MD stuff");
2008 STYPE ("acc-grps", accgrps, NULL);
2009 STYPE ("accelerate", acc, NULL);
2010 STYPE ("freezegrps", freeze, NULL);
2011 STYPE ("freezedim", frdim, NULL);
2012 RTYPE ("cos-acceleration", ir->cos_accel, 0);
2013 STYPE ("deform", deform, NULL);
2015 /* simulated tempering variables */
2016 CCTYPE("simulated tempering variables");
2017 EETYPE("simulated-tempering", ir->bSimTemp, yesno_names);
2018 EETYPE("simulated-tempering-scaling", ir->simtempvals->eSimTempScale, esimtemp_names);
2019 RTYPE("sim-temp-low", ir->simtempvals->simtemp_low, 300.0);
2020 RTYPE("sim-temp-high", ir->simtempvals->simtemp_high, 300.0);
2022 /* expanded ensemble variables */
2023 if (ir->efep == efepEXPANDED || ir->bSimTemp)
2025 read_expandedparams(&ninp, &inp, expand, wi);
2028 /* Electric fields */
2029 CCTYPE("Electric fields");
2030 CTYPE ("Format is number of terms (int) and for all terms an amplitude (real)");
2031 CTYPE ("and a phase angle (real)");
2032 STYPE ("E-x", efield_x, NULL);
2033 STYPE ("E-xt", efield_xt, NULL);
2034 STYPE ("E-y", efield_y, NULL);
2035 STYPE ("E-yt", efield_yt, NULL);
2036 STYPE ("E-z", efield_z, NULL);
2037 STYPE ("E-zt", efield_zt, NULL);
2039 /* AdResS defined thingies */
2040 CCTYPE ("AdResS parameters");
2041 EETYPE("adress", ir->bAdress, yesno_names);
2044 snew(ir->adress, 1);
2045 read_adressparams(&ninp, &inp, ir->adress, wi);
2048 /* User defined thingies */
2049 CCTYPE ("User defined thingies");
2050 STYPE ("user1-grps", user1, NULL);
2051 STYPE ("user2-grps", user2, NULL);
2052 ITYPE ("userint1", ir->userint1, 0);
2053 ITYPE ("userint2", ir->userint2, 0);
2054 ITYPE ("userint3", ir->userint3, 0);
2055 ITYPE ("userint4", ir->userint4, 0);
2056 RTYPE ("userreal1", ir->userreal1, 0);
2057 RTYPE ("userreal2", ir->userreal2, 0);
2058 RTYPE ("userreal3", ir->userreal3, 0);
2059 RTYPE ("userreal4", ir->userreal4, 0);
2062 write_inpfile(mdparout, ninp, inp, FALSE, wi);
2063 for (i = 0; (i < ninp); i++)
2066 sfree(inp[i].value);
2070 /* Process options if necessary */
2071 for (m = 0; m < 2; m++)
2073 for (i = 0; i < 2*DIM; i++)
2082 if (sscanf(dumstr[m], "%lf", &(dumdub[m][XX])) != 1)
2084 warning_error(wi, "Pressure coupling not enough values (I need 1)");
2086 dumdub[m][YY] = dumdub[m][ZZ] = dumdub[m][XX];
2088 case epctSEMIISOTROPIC:
2089 case epctSURFACETENSION:
2090 if (sscanf(dumstr[m], "%lf%lf",
2091 &(dumdub[m][XX]), &(dumdub[m][ZZ])) != 2)
2093 warning_error(wi, "Pressure coupling not enough values (I need 2)");
2095 dumdub[m][YY] = dumdub[m][XX];
2097 case epctANISOTROPIC:
2098 if (sscanf(dumstr[m], "%lf%lf%lf%lf%lf%lf",
2099 &(dumdub[m][XX]), &(dumdub[m][YY]), &(dumdub[m][ZZ]),
2100 &(dumdub[m][3]), &(dumdub[m][4]), &(dumdub[m][5])) != 6)
2102 warning_error(wi, "Pressure coupling not enough values (I need 6)");
2106 gmx_fatal(FARGS, "Pressure coupling type %s not implemented yet",
2107 epcoupltype_names[ir->epct]);
2111 clear_mat(ir->ref_p);
2112 clear_mat(ir->compress);
2113 for (i = 0; i < DIM; i++)
2115 ir->ref_p[i][i] = dumdub[1][i];
2116 ir->compress[i][i] = dumdub[0][i];
2118 if (ir->epct == epctANISOTROPIC)
2120 ir->ref_p[XX][YY] = dumdub[1][3];
2121 ir->ref_p[XX][ZZ] = dumdub[1][4];
2122 ir->ref_p[YY][ZZ] = dumdub[1][5];
2123 if (ir->ref_p[XX][YY] != 0 && ir->ref_p[XX][ZZ] != 0 && ir->ref_p[YY][ZZ] != 0)
2125 warning(wi, "All off-diagonal reference pressures are non-zero. Are you sure you want to apply a threefold shear stress?\n");
2127 ir->compress[XX][YY] = dumdub[0][3];
2128 ir->compress[XX][ZZ] = dumdub[0][4];
2129 ir->compress[YY][ZZ] = dumdub[0][5];
2130 for (i = 0; i < DIM; i++)
2132 for (m = 0; m < i; m++)
2134 ir->ref_p[i][m] = ir->ref_p[m][i];
2135 ir->compress[i][m] = ir->compress[m][i];
2140 if (ir->comm_mode == ecmNO)
2145 opts->couple_moltype = NULL;
2146 if (strlen(couple_moltype) > 0)
2148 if (ir->efep != efepNO)
2150 opts->couple_moltype = strdup(couple_moltype);
2151 if (opts->couple_lam0 == opts->couple_lam1)
2153 warning(wi, "The lambda=0 and lambda=1 states for coupling are identical");
2155 if (ir->eI == eiMD && (opts->couple_lam0 == ecouplamNONE ||
2156 opts->couple_lam1 == ecouplamNONE))
2158 warning(wi, "For proper sampling of the (nearly) decoupled state, stochastic dynamics should be used");
2163 warning(wi, "Can not couple a molecule with free_energy = no");
2166 /* FREE ENERGY AND EXPANDED ENSEMBLE OPTIONS */
2167 if (ir->efep != efepNO)
2169 if (fep->delta_lambda > 0)
2171 ir->efep = efepSLOWGROWTH;
2177 fep->bPrintEnergy = TRUE;
2178 /* always print out the energy to dhdl if we are doing expanded ensemble, since we need the total energy
2179 if the temperature is changing. */
2182 if ((ir->efep != efepNO) || ir->bSimTemp)
2184 ir->bExpanded = FALSE;
2185 if ((ir->efep == efepEXPANDED) || ir->bSimTemp)
2187 ir->bExpanded = TRUE;
2189 do_fep_params(ir, fep_lambda, lambda_weights);
2190 if (ir->bSimTemp) /* done after fep params */
2192 do_simtemp_params(ir);
2197 ir->fepvals->n_lambda = 0;
2200 /* WALL PARAMETERS */
2202 do_wall_params(ir, wall_atomtype, wall_density, opts);
2204 /* ORIENTATION RESTRAINT PARAMETERS */
2206 if (opts->bOrire && str_nelem(orirefitgrp, MAXPTR, NULL) != 1)
2208 warning_error(wi, "ERROR: Need one orientation restraint fit group\n");
2211 /* DEFORMATION PARAMETERS */
2213 clear_mat(ir->deform);
2214 for (i = 0; i < 6; i++)
2218 m = sscanf(deform, "%lf %lf %lf %lf %lf %lf",
2219 &(dumdub[0][0]), &(dumdub[0][1]), &(dumdub[0][2]),
2220 &(dumdub[0][3]), &(dumdub[0][4]), &(dumdub[0][5]));
2221 for (i = 0; i < 3; i++)
2223 ir->deform[i][i] = dumdub[0][i];
2225 ir->deform[YY][XX] = dumdub[0][3];
2226 ir->deform[ZZ][XX] = dumdub[0][4];
2227 ir->deform[ZZ][YY] = dumdub[0][5];
2228 if (ir->epc != epcNO)
2230 for (i = 0; i < 3; i++)
2232 for (j = 0; j <= i; j++)
2234 if (ir->deform[i][j] != 0 && ir->compress[i][j] != 0)
2236 warning_error(wi, "A box element has deform set and compressibility > 0");
2240 for (i = 0; i < 3; i++)
2242 for (j = 0; j < i; j++)
2244 if (ir->deform[i][j] != 0)
2246 for (m = j; m < DIM; m++)
2248 if (ir->compress[m][j] != 0)
2250 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.");
2251 warning(wi, warn_buf);
2263 static int search_QMstring(char *s, int ng, const char *gn[])
2265 /* same as normal search_string, but this one searches QM strings */
2268 for (i = 0; (i < ng); i++)
2270 if (gmx_strcasecmp(s, gn[i]) == 0)
2276 gmx_fatal(FARGS, "this QM method or basisset (%s) is not implemented\n!", s);
2280 } /* search_QMstring */
2283 int search_string(char *s, int ng, char *gn[])
2287 for (i = 0; (i < ng); i++)
2289 if (gmx_strcasecmp(s, gn[i]) == 0)
2296 "Group %s referenced in the .mdp file was not found in the index file.\n"
2297 "Group names must match either [moleculetype] names or custom index group\n"
2298 "names, in which case you must supply an index file to the '-n' option\n"
2305 static gmx_bool do_numbering(int natoms, gmx_groups_t *groups, int ng, char *ptrs[],
2306 t_blocka *block, char *gnames[],
2307 int gtype, int restnm,
2308 int grptp, gmx_bool bVerbose,
2311 unsigned short *cbuf;
2312 t_grps *grps = &(groups->grps[gtype]);
2313 int i, j, gid, aj, ognr, ntot = 0;
2316 char warn_buf[STRLEN];
2320 fprintf(debug, "Starting numbering %d groups of type %d\n", ng, gtype);
2323 title = gtypes[gtype];
2326 /* Mark all id's as not set */
2327 for (i = 0; (i < natoms); i++)
2332 snew(grps->nm_ind, ng+1); /* +1 for possible rest group */
2333 for (i = 0; (i < ng); i++)
2335 /* Lookup the group name in the block structure */
2336 gid = search_string(ptrs[i], block->nr, gnames);
2337 if ((grptp != egrptpONE) || (i == 0))
2339 grps->nm_ind[grps->nr++] = gid;
2343 fprintf(debug, "Found gid %d for group %s\n", gid, ptrs[i]);
2346 /* Now go over the atoms in the group */
2347 for (j = block->index[gid]; (j < block->index[gid+1]); j++)
2352 /* Range checking */
2353 if ((aj < 0) || (aj >= natoms))
2355 gmx_fatal(FARGS, "Invalid atom number %d in indexfile", aj);
2357 /* Lookup up the old group number */
2361 gmx_fatal(FARGS, "Atom %d in multiple %s groups (%d and %d)",
2362 aj+1, title, ognr+1, i+1);
2366 /* Store the group number in buffer */
2367 if (grptp == egrptpONE)
2380 /* Now check whether we have done all atoms */
2384 if (grptp == egrptpALL)
2386 gmx_fatal(FARGS, "%d atoms are not part of any of the %s groups",
2387 natoms-ntot, title);
2389 else if (grptp == egrptpPART)
2391 sprintf(warn_buf, "%d atoms are not part of any of the %s groups",
2392 natoms-ntot, title);
2393 warning_note(wi, warn_buf);
2395 /* Assign all atoms currently unassigned to a rest group */
2396 for (j = 0; (j < natoms); j++)
2398 if (cbuf[j] == NOGID)
2404 if (grptp != egrptpPART)
2409 "Making dummy/rest group for %s containing %d elements\n",
2410 title, natoms-ntot);
2412 /* Add group name "rest" */
2413 grps->nm_ind[grps->nr] = restnm;
2415 /* Assign the rest name to all atoms not currently assigned to a group */
2416 for (j = 0; (j < natoms); j++)
2418 if (cbuf[j] == NOGID)
2427 if (grps->nr == 1 && (ntot == 0 || ntot == natoms))
2429 /* All atoms are part of one (or no) group, no index required */
2430 groups->ngrpnr[gtype] = 0;
2431 groups->grpnr[gtype] = NULL;
2435 groups->ngrpnr[gtype] = natoms;
2436 snew(groups->grpnr[gtype], natoms);
2437 for (j = 0; (j < natoms); j++)
2439 groups->grpnr[gtype][j] = cbuf[j];
2445 return (bRest && grptp == egrptpPART);
2448 static void calc_nrdf(gmx_mtop_t *mtop, t_inputrec *ir, char **gnames)
2451 gmx_groups_t *groups;
2453 int natoms, ai, aj, i, j, d, g, imin, jmin, nc;
2455 int *nrdf2, *na_vcm, na_tot;
2456 double *nrdf_tc, *nrdf_vcm, nrdf_uc, n_sub = 0;
2457 gmx_mtop_atomloop_all_t aloop;
2459 int mb, mol, ftype, as;
2460 gmx_molblock_t *molb;
2461 gmx_moltype_t *molt;
2464 * First calc 3xnr-atoms for each group
2465 * then subtract half a degree of freedom for each constraint
2467 * Only atoms and nuclei contribute to the degrees of freedom...
2472 groups = &mtop->groups;
2473 natoms = mtop->natoms;
2475 /* Allocate one more for a possible rest group */
2476 /* We need to sum degrees of freedom into doubles,
2477 * since floats give too low nrdf's above 3 million atoms.
2479 snew(nrdf_tc, groups->grps[egcTC].nr+1);
2480 snew(nrdf_vcm, groups->grps[egcVCM].nr+1);
2481 snew(na_vcm, groups->grps[egcVCM].nr+1);
2483 for (i = 0; i < groups->grps[egcTC].nr; i++)
2487 for (i = 0; i < groups->grps[egcVCM].nr+1; i++)
2492 snew(nrdf2, natoms);
2493 aloop = gmx_mtop_atomloop_all_init(mtop);
2494 while (gmx_mtop_atomloop_all_next(aloop, &i, &atom))
2497 if (atom->ptype == eptAtom || atom->ptype == eptNucleus)
2499 g = ggrpnr(groups, egcFREEZE, i);
2500 /* Double count nrdf for particle i */
2501 for (d = 0; d < DIM; d++)
2503 if (opts->nFreeze[g][d] == 0)
2508 nrdf_tc [ggrpnr(groups, egcTC, i)] += 0.5*nrdf2[i];
2509 nrdf_vcm[ggrpnr(groups, egcVCM, i)] += 0.5*nrdf2[i];
2514 for (mb = 0; mb < mtop->nmolblock; mb++)
2516 molb = &mtop->molblock[mb];
2517 molt = &mtop->moltype[molb->type];
2518 atom = molt->atoms.atom;
2519 for (mol = 0; mol < molb->nmol; mol++)
2521 for (ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++)
2523 ia = molt->ilist[ftype].iatoms;
2524 for (i = 0; i < molt->ilist[ftype].nr; )
2526 /* Subtract degrees of freedom for the constraints,
2527 * if the particles still have degrees of freedom left.
2528 * If one of the particles is a vsite or a shell, then all
2529 * constraint motion will go there, but since they do not
2530 * contribute to the constraints the degrees of freedom do not
2535 if (((atom[ia[1]].ptype == eptNucleus) ||
2536 (atom[ia[1]].ptype == eptAtom)) &&
2537 ((atom[ia[2]].ptype == eptNucleus) ||
2538 (atom[ia[2]].ptype == eptAtom)))
2556 imin = min(imin, nrdf2[ai]);
2557 jmin = min(jmin, nrdf2[aj]);
2560 nrdf_tc [ggrpnr(groups, egcTC, ai)] -= 0.5*imin;
2561 nrdf_tc [ggrpnr(groups, egcTC, aj)] -= 0.5*jmin;
2562 nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5*imin;
2563 nrdf_vcm[ggrpnr(groups, egcVCM, aj)] -= 0.5*jmin;
2565 ia += interaction_function[ftype].nratoms+1;
2566 i += interaction_function[ftype].nratoms+1;
2569 ia = molt->ilist[F_SETTLE].iatoms;
2570 for (i = 0; i < molt->ilist[F_SETTLE].nr; )
2572 /* Subtract 1 dof from every atom in the SETTLE */
2573 for (j = 0; j < 3; j++)
2576 imin = min(2, nrdf2[ai]);
2578 nrdf_tc [ggrpnr(groups, egcTC, ai)] -= 0.5*imin;
2579 nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5*imin;
2584 as += molt->atoms.nr;
2588 if (ir->ePull == epullCONSTRAINT)
2590 /* Correct nrdf for the COM constraints.
2591 * We correct using the TC and VCM group of the first atom
2592 * in the reference and pull group. If atoms in one pull group
2593 * belong to different TC or VCM groups it is anyhow difficult
2594 * to determine the optimal nrdf assignment.
2597 if (pull->eGeom == epullgPOS)
2600 for (i = 0; i < DIM; i++)
2612 for (i = 0; i < pull->ngrp; i++)
2615 if (pull->grp[0].nat > 0)
2617 /* Subtract 1/2 dof from the reference group */
2618 ai = pull->grp[0].ind[0];
2619 if (nrdf_tc[ggrpnr(groups, egcTC, ai)] > 1)
2621 nrdf_tc [ggrpnr(groups, egcTC, ai)] -= 0.5;
2622 nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5;
2626 /* Subtract 1/2 dof from the pulled group */
2627 ai = pull->grp[1+i].ind[0];
2628 nrdf_tc [ggrpnr(groups, egcTC, ai)] -= 0.5*imin;
2629 nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5*imin;
2630 if (nrdf_tc[ggrpnr(groups, egcTC, ai)] < 0)
2632 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)]]);
2637 if (ir->nstcomm != 0)
2639 /* Subtract 3 from the number of degrees of freedom in each vcm group
2640 * when com translation is removed and 6 when rotation is removed
2643 switch (ir->comm_mode)
2646 n_sub = ndof_com(ir);
2653 gmx_incons("Checking comm_mode");
2656 for (i = 0; i < groups->grps[egcTC].nr; i++)
2658 /* Count the number of atoms of TC group i for every VCM group */
2659 for (j = 0; j < groups->grps[egcVCM].nr+1; j++)
2664 for (ai = 0; ai < natoms; ai++)
2666 if (ggrpnr(groups, egcTC, ai) == i)
2668 na_vcm[ggrpnr(groups, egcVCM, ai)]++;
2672 /* Correct for VCM removal according to the fraction of each VCM
2673 * group present in this TC group.
2675 nrdf_uc = nrdf_tc[i];
2678 fprintf(debug, "T-group[%d] nrdf_uc = %g, n_sub = %g\n",
2682 for (j = 0; j < groups->grps[egcVCM].nr+1; j++)
2684 if (nrdf_vcm[j] > n_sub)
2686 nrdf_tc[i] += nrdf_uc*((double)na_vcm[j]/(double)na_tot)*
2687 (nrdf_vcm[j] - n_sub)/nrdf_vcm[j];
2691 fprintf(debug, " nrdf_vcm[%d] = %g, nrdf = %g\n",
2692 j, nrdf_vcm[j], nrdf_tc[i]);
2697 for (i = 0; (i < groups->grps[egcTC].nr); i++)
2699 opts->nrdf[i] = nrdf_tc[i];
2700 if (opts->nrdf[i] < 0)
2705 "Number of degrees of freedom in T-Coupling group %s is %.2f\n",
2706 gnames[groups->grps[egcTC].nm_ind[i]], opts->nrdf[i]);
2715 static void decode_cos(char *s, t_cosines *cosine, gmx_bool bTime)
2718 char format[STRLEN], f1[STRLEN];
2730 sscanf(t, "%d", &(cosine->n));
2737 snew(cosine->a, cosine->n);
2738 snew(cosine->phi, cosine->n);
2740 sprintf(format, "%%*d");
2741 for (i = 0; (i < cosine->n); i++)
2744 strcat(f1, "%lf%lf");
2745 if (sscanf(t, f1, &a, &phi) < 2)
2747 gmx_fatal(FARGS, "Invalid input for electric field shift: '%s'", t);
2750 cosine->phi[i] = phi;
2751 strcat(format, "%*lf%*lf");
2758 static gmx_bool do_egp_flag(t_inputrec *ir, gmx_groups_t *groups,
2759 const char *option, const char *val, int flag)
2761 /* The maximum number of energy group pairs would be MAXPTR*(MAXPTR+1)/2.
2762 * But since this is much larger than STRLEN, such a line can not be parsed.
2763 * The real maximum is the number of names that fit in a string: STRLEN/2.
2765 #define EGP_MAX (STRLEN/2)
2766 int nelem, i, j, k, nr;
2767 char *names[EGP_MAX];
2771 gnames = groups->grpname;
2773 nelem = str_nelem(val, EGP_MAX, names);
2776 gmx_fatal(FARGS, "The number of groups for %s is odd", option);
2778 nr = groups->grps[egcENER].nr;
2780 for (i = 0; i < nelem/2; i++)
2784 gmx_strcasecmp(names[2*i], *(gnames[groups->grps[egcENER].nm_ind[j]])))
2790 gmx_fatal(FARGS, "%s in %s is not an energy group\n",
2791 names[2*i], option);
2795 gmx_strcasecmp(names[2*i+1], *(gnames[groups->grps[egcENER].nm_ind[k]])))
2801 gmx_fatal(FARGS, "%s in %s is not an energy group\n",
2802 names[2*i+1], option);
2804 if ((j < nr) && (k < nr))
2806 ir->opts.egp_flags[nr*j+k] |= flag;
2807 ir->opts.egp_flags[nr*k+j] |= flag;
2815 void do_index(const char* mdparin, const char *ndx,
2818 t_inputrec *ir, rvec *v,
2822 gmx_groups_t *groups;
2826 char warnbuf[STRLEN], **gnames;
2827 int nr, ntcg, ntau_t, nref_t, nacc, nofg, nSA, nSA_points, nSA_time, nSA_temp;
2830 int nacg, nfreeze, nfrdim, nenergy, nvcm, nuser;
2831 char *ptr1[MAXPTR], *ptr2[MAXPTR], *ptr3[MAXPTR];
2832 int i, j, k, restnm;
2834 gmx_bool bExcl, bTable, bSetTCpar, bAnneal, bRest;
2835 int nQMmethod, nQMbasis, nQMcharge, nQMmult, nbSH, nCASorb, nCASelec,
2836 nSAon, nSAoff, nSAsteps, nQMg, nbOPT, nbTS;
2837 char warn_buf[STRLEN];
2841 fprintf(stderr, "processing index file...\n");
2847 snew(grps->index, 1);
2849 atoms_all = gmx_mtop_global_atoms(mtop);
2850 analyse(&atoms_all, grps, &gnames, FALSE, TRUE);
2851 free_t_atoms(&atoms_all, FALSE);
2855 grps = init_index(ndx, &gnames);
2858 groups = &mtop->groups;
2859 natoms = mtop->natoms;
2860 symtab = &mtop->symtab;
2862 snew(groups->grpname, grps->nr+1);
2864 for (i = 0; (i < grps->nr); i++)
2866 groups->grpname[i] = put_symtab(symtab, gnames[i]);
2868 groups->grpname[i] = put_symtab(symtab, "rest");
2870 srenew(gnames, grps->nr+1);
2871 gnames[restnm] = *(groups->grpname[i]);
2872 groups->ngrpname = grps->nr+1;
2874 set_warning_line(wi, mdparin, -1);
2876 ntau_t = str_nelem(tau_t, MAXPTR, ptr1);
2877 nref_t = str_nelem(ref_t, MAXPTR, ptr2);
2878 ntcg = str_nelem(tcgrps, MAXPTR, ptr3);
2879 if ((ntau_t != ntcg) || (nref_t != ntcg))
2881 gmx_fatal(FARGS, "Invalid T coupling input: %d groups, %d ref-t values and "
2882 "%d tau-t values", ntcg, nref_t, ntau_t);
2885 bSetTCpar = (ir->etc || EI_SD(ir->eI) || ir->eI == eiBD || EI_TPI(ir->eI));
2886 do_numbering(natoms, groups, ntcg, ptr3, grps, gnames, egcTC,
2887 restnm, bSetTCpar ? egrptpALL : egrptpALL_GENREST, bVerbose, wi);
2888 nr = groups->grps[egcTC].nr;
2890 snew(ir->opts.nrdf, nr);
2891 snew(ir->opts.tau_t, nr);
2892 snew(ir->opts.ref_t, nr);
2893 if (ir->eI == eiBD && ir->bd_fric == 0)
2895 fprintf(stderr, "bd-fric=0, so tau-t will be used as the inverse friction constant(s)\n");
2902 gmx_fatal(FARGS, "Not enough ref-t and tau-t values!");
2906 for (i = 0; (i < nr); i++)
2908 ir->opts.tau_t[i] = strtod(ptr1[i], NULL);
2909 if ((ir->eI == eiBD || ir->eI == eiSD2) && ir->opts.tau_t[i] <= 0)
2911 sprintf(warn_buf, "With integrator %s tau-t should be larger than 0", ei_names[ir->eI]);
2912 warning_error(wi, warn_buf);
2915 if (ir->etc != etcVRESCALE && ir->opts.tau_t[i] == 0)
2917 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.");
2920 if (ir->opts.tau_t[i] >= 0)
2922 tau_min = min(tau_min, ir->opts.tau_t[i]);
2925 if (ir->etc != etcNO && ir->nsttcouple == -1)
2927 ir->nsttcouple = ir_optimal_nsttcouple(ir);
2932 if ((ir->etc == etcNOSEHOOVER) && (ir->epc == epcBERENDSEN))
2934 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");
2936 if ((ir->epc == epcMTTK) && (ir->etc > etcNO))
2938 if (ir->nstpcouple != ir->nsttcouple)
2940 int mincouple = min(ir->nstpcouple, ir->nsttcouple);
2941 ir->nstpcouple = ir->nsttcouple = mincouple;
2942 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);
2943 warning_note(wi, warn_buf);
2947 /* velocity verlet with averaged kinetic energy KE = 0.5*(v(t+1/2) - v(t-1/2)) is implemented
2948 primarily for testing purposes, and does not work with temperature coupling other than 1 */
2950 if (ETC_ANDERSEN(ir->etc))
2952 if (ir->nsttcouple != 1)
2955 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");
2956 warning_note(wi, warn_buf);
2959 nstcmin = tcouple_min_integration_steps(ir->etc);
2962 if (tau_min/(ir->delta_t*ir->nsttcouple) < nstcmin)
2964 sprintf(warn_buf, "For proper integration of the %s thermostat, tau-t (%g) should be at least %d times larger than nsttcouple*dt (%g)",
2965 ETCOUPLTYPE(ir->etc),
2967 ir->nsttcouple*ir->delta_t);
2968 warning(wi, warn_buf);
2971 for (i = 0; (i < nr); i++)
2973 ir->opts.ref_t[i] = strtod(ptr2[i], NULL);
2974 if (ir->opts.ref_t[i] < 0)
2976 gmx_fatal(FARGS, "ref-t for group %d negative", i);
2979 /* set the lambda mc temperature to the md integrator temperature (which should be defined
2980 if we are in this conditional) if mc_temp is negative */
2981 if (ir->expandedvals->mc_temp < 0)
2983 ir->expandedvals->mc_temp = ir->opts.ref_t[0]; /*for now, set to the first reft */
2987 /* Simulated annealing for each group. There are nr groups */
2988 nSA = str_nelem(anneal, MAXPTR, ptr1);
2989 if (nSA == 1 && (ptr1[0][0] == 'n' || ptr1[0][0] == 'N'))
2993 if (nSA > 0 && nSA != nr)
2995 gmx_fatal(FARGS, "Not enough annealing values: %d (for %d groups)\n", nSA, nr);
2999 snew(ir->opts.annealing, nr);
3000 snew(ir->opts.anneal_npoints, nr);
3001 snew(ir->opts.anneal_time, nr);
3002 snew(ir->opts.anneal_temp, nr);
3003 for (i = 0; i < nr; i++)
3005 ir->opts.annealing[i] = eannNO;
3006 ir->opts.anneal_npoints[i] = 0;
3007 ir->opts.anneal_time[i] = NULL;
3008 ir->opts.anneal_temp[i] = NULL;
3013 for (i = 0; i < nr; i++)
3015 if (ptr1[i][0] == 'n' || ptr1[i][0] == 'N')
3017 ir->opts.annealing[i] = eannNO;
3019 else if (ptr1[i][0] == 's' || ptr1[i][0] == 'S')
3021 ir->opts.annealing[i] = eannSINGLE;
3024 else if (ptr1[i][0] == 'p' || ptr1[i][0] == 'P')
3026 ir->opts.annealing[i] = eannPERIODIC;
3032 /* Read the other fields too */
3033 nSA_points = str_nelem(anneal_npoints, MAXPTR, ptr1);
3034 if (nSA_points != nSA)
3036 gmx_fatal(FARGS, "Found %d annealing-npoints values for %d groups\n", nSA_points, nSA);
3038 for (k = 0, i = 0; i < nr; i++)
3040 ir->opts.anneal_npoints[i] = strtol(ptr1[i], NULL, 10);
3041 if (ir->opts.anneal_npoints[i] == 1)
3043 gmx_fatal(FARGS, "Please specify at least a start and an end point for annealing\n");
3045 snew(ir->opts.anneal_time[i], ir->opts.anneal_npoints[i]);
3046 snew(ir->opts.anneal_temp[i], ir->opts.anneal_npoints[i]);
3047 k += ir->opts.anneal_npoints[i];
3050 nSA_time = str_nelem(anneal_time, MAXPTR, ptr1);
3053 gmx_fatal(FARGS, "Found %d annealing-time values, wanter %d\n", nSA_time, k);
3055 nSA_temp = str_nelem(anneal_temp, MAXPTR, ptr2);
3058 gmx_fatal(FARGS, "Found %d annealing-temp values, wanted %d\n", nSA_temp, k);
3061 for (i = 0, k = 0; i < nr; i++)
3064 for (j = 0; j < ir->opts.anneal_npoints[i]; j++)
3066 ir->opts.anneal_time[i][j] = strtod(ptr1[k], NULL);
3067 ir->opts.anneal_temp[i][j] = strtod(ptr2[k], NULL);
3070 if (ir->opts.anneal_time[i][0] > (ir->init_t+GMX_REAL_EPS))
3072 gmx_fatal(FARGS, "First time point for annealing > init_t.\n");
3078 if (ir->opts.anneal_time[i][j] < ir->opts.anneal_time[i][j-1])
3080 gmx_fatal(FARGS, "Annealing timepoints out of order: t=%f comes after t=%f\n",
3081 ir->opts.anneal_time[i][j], ir->opts.anneal_time[i][j-1]);
3084 if (ir->opts.anneal_temp[i][j] < 0)
3086 gmx_fatal(FARGS, "Found negative temperature in annealing: %f\n", ir->opts.anneal_temp[i][j]);
3091 /* Print out some summary information, to make sure we got it right */
3092 for (i = 0, k = 0; i < nr; i++)
3094 if (ir->opts.annealing[i] != eannNO)
3096 j = groups->grps[egcTC].nm_ind[i];
3097 fprintf(stderr, "Simulated annealing for group %s: %s, %d timepoints\n",
3098 *(groups->grpname[j]), eann_names[ir->opts.annealing[i]],
3099 ir->opts.anneal_npoints[i]);
3100 fprintf(stderr, "Time (ps) Temperature (K)\n");
3101 /* All terms except the last one */
3102 for (j = 0; j < (ir->opts.anneal_npoints[i]-1); j++)
3104 fprintf(stderr, "%9.1f %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3107 /* Finally the last one */
3108 j = ir->opts.anneal_npoints[i]-1;
3109 if (ir->opts.annealing[i] == eannSINGLE)
3111 fprintf(stderr, "%9.1f- %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3115 fprintf(stderr, "%9.1f %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3116 if (fabs(ir->opts.anneal_temp[i][j]-ir->opts.anneal_temp[i][0]) > GMX_REAL_EPS)
3118 warning_note(wi, "There is a temperature jump when your annealing loops back.\n");
3127 if (ir->ePull != epullNO)
3129 make_pull_groups(ir->pull, pull_grp, grps, gnames);
3134 make_rotation_groups(ir->rot, rot_grp, grps, gnames);
3137 nacc = str_nelem(acc, MAXPTR, ptr1);
3138 nacg = str_nelem(accgrps, MAXPTR, ptr2);
3139 if (nacg*DIM != nacc)
3141 gmx_fatal(FARGS, "Invalid Acceleration input: %d groups and %d acc. values",
3144 do_numbering(natoms, groups, nacg, ptr2, grps, gnames, egcACC,
3145 restnm, egrptpALL_GENREST, bVerbose, wi);
3146 nr = groups->grps[egcACC].nr;
3147 snew(ir->opts.acc, nr);
3148 ir->opts.ngacc = nr;
3150 for (i = k = 0; (i < nacg); i++)
3152 for (j = 0; (j < DIM); j++, k++)
3154 ir->opts.acc[i][j] = strtod(ptr1[k], NULL);
3157 for (; (i < nr); i++)
3159 for (j = 0; (j < DIM); j++)
3161 ir->opts.acc[i][j] = 0;
3165 nfrdim = str_nelem(frdim, MAXPTR, ptr1);
3166 nfreeze = str_nelem(freeze, MAXPTR, ptr2);
3167 if (nfrdim != DIM*nfreeze)
3169 gmx_fatal(FARGS, "Invalid Freezing input: %d groups and %d freeze values",
3172 do_numbering(natoms, groups, nfreeze, ptr2, grps, gnames, egcFREEZE,
3173 restnm, egrptpALL_GENREST, bVerbose, wi);
3174 nr = groups->grps[egcFREEZE].nr;
3175 ir->opts.ngfrz = nr;
3176 snew(ir->opts.nFreeze, nr);
3177 for (i = k = 0; (i < nfreeze); i++)
3179 for (j = 0; (j < DIM); j++, k++)
3181 ir->opts.nFreeze[i][j] = (gmx_strncasecmp(ptr1[k], "Y", 1) == 0);
3182 if (!ir->opts.nFreeze[i][j])
3184 if (gmx_strncasecmp(ptr1[k], "N", 1) != 0)
3186 sprintf(warnbuf, "Please use Y(ES) or N(O) for freezedim only "
3187 "(not %s)", ptr1[k]);
3188 warning(wi, warn_buf);
3193 for (; (i < nr); i++)
3195 for (j = 0; (j < DIM); j++)
3197 ir->opts.nFreeze[i][j] = 0;
3201 nenergy = str_nelem(energy, MAXPTR, ptr1);
3202 do_numbering(natoms, groups, nenergy, ptr1, grps, gnames, egcENER,
3203 restnm, egrptpALL_GENREST, bVerbose, wi);
3204 add_wall_energrps(groups, ir->nwall, symtab);
3205 ir->opts.ngener = groups->grps[egcENER].nr;
3206 nvcm = str_nelem(vcm, MAXPTR, ptr1);
3208 do_numbering(natoms, groups, nvcm, ptr1, grps, gnames, egcVCM,
3209 restnm, nvcm == 0 ? egrptpALL_GENREST : egrptpPART, bVerbose, wi);
3212 warning(wi, "Some atoms are not part of any center of mass motion removal group.\n"
3213 "This may lead to artifacts.\n"
3214 "In most cases one should use one group for the whole system.");
3217 /* Now we have filled the freeze struct, so we can calculate NRDF */
3218 calc_nrdf(mtop, ir, gnames);
3224 /* Must check per group! */
3225 for (i = 0; (i < ir->opts.ngtc); i++)
3227 ntot += ir->opts.nrdf[i];
3229 if (ntot != (DIM*natoms))
3231 fac = sqrt(ntot/(DIM*natoms));
3234 fprintf(stderr, "Scaling velocities by a factor of %.3f to account for constraints\n"
3235 "and removal of center of mass motion\n", fac);
3237 for (i = 0; (i < natoms); i++)
3239 svmul(fac, v[i], v[i]);
3244 nuser = str_nelem(user1, MAXPTR, ptr1);
3245 do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcUser1,
3246 restnm, egrptpALL_GENREST, bVerbose, wi);
3247 nuser = str_nelem(user2, MAXPTR, ptr1);
3248 do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcUser2,
3249 restnm, egrptpALL_GENREST, bVerbose, wi);
3250 nuser = str_nelem(xtc_grps, MAXPTR, ptr1);
3251 do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcXTC,
3252 restnm, egrptpONE, bVerbose, wi);
3253 nofg = str_nelem(orirefitgrp, MAXPTR, ptr1);
3254 do_numbering(natoms, groups, nofg, ptr1, grps, gnames, egcORFIT,
3255 restnm, egrptpALL_GENREST, bVerbose, wi);
3257 /* QMMM input processing */
3258 nQMg = str_nelem(QMMM, MAXPTR, ptr1);
3259 nQMmethod = str_nelem(QMmethod, MAXPTR, ptr2);
3260 nQMbasis = str_nelem(QMbasis, MAXPTR, ptr3);
3261 if ((nQMmethod != nQMg) || (nQMbasis != nQMg))
3263 gmx_fatal(FARGS, "Invalid QMMM input: %d groups %d basissets"
3264 " and %d methods\n", nQMg, nQMbasis, nQMmethod);
3266 /* group rest, if any, is always MM! */
3267 do_numbering(natoms, groups, nQMg, ptr1, grps, gnames, egcQMMM,
3268 restnm, egrptpALL_GENREST, bVerbose, wi);
3269 nr = nQMg; /*atoms->grps[egcQMMM].nr;*/
3270 ir->opts.ngQM = nQMg;
3271 snew(ir->opts.QMmethod, nr);
3272 snew(ir->opts.QMbasis, nr);
3273 for (i = 0; i < nr; i++)
3275 /* input consists of strings: RHF CASSCF PM3 .. These need to be
3276 * converted to the corresponding enum in names.c
3278 ir->opts.QMmethod[i] = search_QMstring(ptr2[i], eQMmethodNR,
3280 ir->opts.QMbasis[i] = search_QMstring(ptr3[i], eQMbasisNR,
3284 nQMmult = str_nelem(QMmult, MAXPTR, ptr1);
3285 nQMcharge = str_nelem(QMcharge, MAXPTR, ptr2);
3286 nbSH = str_nelem(bSH, MAXPTR, ptr3);
3287 snew(ir->opts.QMmult, nr);
3288 snew(ir->opts.QMcharge, nr);
3289 snew(ir->opts.bSH, nr);
3291 for (i = 0; i < nr; i++)
3293 ir->opts.QMmult[i] = strtol(ptr1[i], NULL, 10);
3294 ir->opts.QMcharge[i] = strtol(ptr2[i], NULL, 10);
3295 ir->opts.bSH[i] = (gmx_strncasecmp(ptr3[i], "Y", 1) == 0);
3298 nCASelec = str_nelem(CASelectrons, MAXPTR, ptr1);
3299 nCASorb = str_nelem(CASorbitals, MAXPTR, ptr2);
3300 snew(ir->opts.CASelectrons, nr);
3301 snew(ir->opts.CASorbitals, nr);
3302 for (i = 0; i < nr; i++)
3304 ir->opts.CASelectrons[i] = strtol(ptr1[i], NULL, 10);
3305 ir->opts.CASorbitals[i] = strtol(ptr2[i], NULL, 10);
3307 /* special optimization options */
3309 nbOPT = str_nelem(bOPT, MAXPTR, ptr1);
3310 nbTS = str_nelem(bTS, MAXPTR, ptr2);
3311 snew(ir->opts.bOPT, nr);
3312 snew(ir->opts.bTS, nr);
3313 for (i = 0; i < nr; i++)
3315 ir->opts.bOPT[i] = (gmx_strncasecmp(ptr1[i], "Y", 1) == 0);
3316 ir->opts.bTS[i] = (gmx_strncasecmp(ptr2[i], "Y", 1) == 0);
3318 nSAon = str_nelem(SAon, MAXPTR, ptr1);
3319 nSAoff = str_nelem(SAoff, MAXPTR, ptr2);
3320 nSAsteps = str_nelem(SAsteps, MAXPTR, ptr3);
3321 snew(ir->opts.SAon, nr);
3322 snew(ir->opts.SAoff, nr);
3323 snew(ir->opts.SAsteps, nr);
3325 for (i = 0; i < nr; i++)
3327 ir->opts.SAon[i] = strtod(ptr1[i], NULL);
3328 ir->opts.SAoff[i] = strtod(ptr2[i], NULL);
3329 ir->opts.SAsteps[i] = strtol(ptr3[i], NULL, 10);
3331 /* end of QMMM input */
3335 for (i = 0; (i < egcNR); i++)
3337 fprintf(stderr, "%-16s has %d element(s):", gtypes[i], groups->grps[i].nr);
3338 for (j = 0; (j < groups->grps[i].nr); j++)
3340 fprintf(stderr, " %s", *(groups->grpname[groups->grps[i].nm_ind[j]]));
3342 fprintf(stderr, "\n");
3346 nr = groups->grps[egcENER].nr;
3347 snew(ir->opts.egp_flags, nr*nr);
3349 bExcl = do_egp_flag(ir, groups, "energygrp-excl", egpexcl, EGP_EXCL);
3350 if (bExcl && ir->cutoff_scheme == ecutsVERLET)
3352 warning_error(wi, "Energy group exclusions are not (yet) implemented for the Verlet scheme");
3354 if (bExcl && EEL_FULL(ir->coulombtype))
3356 warning(wi, "Can not exclude the lattice Coulomb energy between energy groups");
3359 bTable = do_egp_flag(ir, groups, "energygrp-table", egptable, EGP_TABLE);
3360 if (bTable && !(ir->vdwtype == evdwUSER) &&
3361 !(ir->coulombtype == eelUSER) && !(ir->coulombtype == eelPMEUSER) &&
3362 !(ir->coulombtype == eelPMEUSERSWITCH))
3364 gmx_fatal(FARGS, "Can only have energy group pair tables in combination with user tables for VdW and/or Coulomb");
3367 decode_cos(efield_x, &(ir->ex[XX]), FALSE);
3368 decode_cos(efield_xt, &(ir->et[XX]), TRUE);
3369 decode_cos(efield_y, &(ir->ex[YY]), FALSE);
3370 decode_cos(efield_yt, &(ir->et[YY]), TRUE);
3371 decode_cos(efield_z, &(ir->ex[ZZ]), FALSE);
3372 decode_cos(efield_zt, &(ir->et[ZZ]), TRUE);
3376 do_adress_index(ir->adress, groups, gnames, &(ir->opts), wi);
3379 for (i = 0; (i < grps->nr); i++)
3391 static void check_disre(gmx_mtop_t *mtop)
3393 gmx_ffparams_t *ffparams;
3394 t_functype *functype;
3396 int i, ndouble, ftype;
3397 int label, old_label;
3399 if (gmx_mtop_ftype_count(mtop, F_DISRES) > 0)
3401 ffparams = &mtop->ffparams;
3402 functype = ffparams->functype;
3403 ip = ffparams->iparams;
3406 for (i = 0; i < ffparams->ntypes; i++)
3408 ftype = functype[i];
3409 if (ftype == F_DISRES)
3411 label = ip[i].disres.label;
3412 if (label == old_label)
3414 fprintf(stderr, "Distance restraint index %d occurs twice\n", label);
3422 gmx_fatal(FARGS, "Found %d double distance restraint indices,\n"
3423 "probably the parameters for multiple pairs in one restraint "
3424 "are not identical\n", ndouble);
3429 static gmx_bool absolute_reference(t_inputrec *ir, gmx_mtop_t *sys,
3430 gmx_bool posres_only,
3434 gmx_mtop_ilistloop_t iloop;
3444 for (d = 0; d < DIM; d++)
3446 AbsRef[d] = (d < ndof_com(ir) ? 0 : 1);
3448 /* Check for freeze groups */
3449 for (g = 0; g < ir->opts.ngfrz; g++)
3451 for (d = 0; d < DIM; d++)
3453 if (ir->opts.nFreeze[g][d] != 0)
3461 /* Check for position restraints */
3462 iloop = gmx_mtop_ilistloop_init(sys);
3463 while (gmx_mtop_ilistloop_next(iloop, &ilist, &nmol))
3466 (AbsRef[XX] == 0 || AbsRef[YY] == 0 || AbsRef[ZZ] == 0))
3468 for (i = 0; i < ilist[F_POSRES].nr; i += 2)
3470 pr = &sys->ffparams.iparams[ilist[F_POSRES].iatoms[i]];
3471 for (d = 0; d < DIM; d++)
3473 if (pr->posres.fcA[d] != 0)
3479 for (i = 0; i < ilist[F_FBPOSRES].nr; i += 2)
3481 /* Check for flat-bottom posres */
3482 pr = &sys->ffparams.iparams[ilist[F_FBPOSRES].iatoms[i]];
3483 if (pr->fbposres.k != 0)
3485 switch (pr->fbposres.geom)
3487 case efbposresSPHERE:
3488 AbsRef[XX] = AbsRef[YY] = AbsRef[ZZ] = 1;
3490 case efbposresCYLINDER:
3491 AbsRef[XX] = AbsRef[YY] = 1;
3493 case efbposresX: /* d=XX */
3494 case efbposresY: /* d=YY */
3495 case efbposresZ: /* d=ZZ */
3496 d = pr->fbposres.geom - efbposresX;
3500 gmx_fatal(FARGS, " Invalid geometry for flat-bottom position restraint.\n"
3501 "Expected nr between 1 and %d. Found %d\n", efbposresNR-1,
3509 return (AbsRef[XX] != 0 && AbsRef[YY] != 0 && AbsRef[ZZ] != 0);
3512 void triple_check(const char *mdparin, t_inputrec *ir, gmx_mtop_t *sys,
3516 int i, m, g, nmol, npct;
3517 gmx_bool bCharge, bAcc;
3518 real gdt_max, *mgrp, mt;
3520 gmx_mtop_atomloop_block_t aloopb;
3521 gmx_mtop_atomloop_all_t aloop;
3524 char warn_buf[STRLEN];
3526 set_warning_line(wi, mdparin, -1);
3528 if (EI_DYNAMICS(ir->eI) && !EI_SD(ir->eI) && ir->eI != eiBD &&
3529 ir->comm_mode == ecmNO &&
3530 !(absolute_reference(ir, sys, FALSE, AbsRef) || ir->nsteps <= 10))
3532 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");
3535 /* Check for pressure coupling with absolute position restraints */
3536 if (ir->epc != epcNO && ir->refcoord_scaling == erscNO)
3538 absolute_reference(ir, sys, TRUE, AbsRef);
3540 for (m = 0; m < DIM; m++)
3542 if (AbsRef[m] && norm2(ir->compress[m]) > 0)
3544 warning(wi, "You are using pressure coupling with absolute position restraints, this will give artifacts. Use the refcoord_scaling option.");
3552 aloopb = gmx_mtop_atomloop_block_init(sys);
3553 while (gmx_mtop_atomloop_block_next(aloopb, &atom, &nmol))
3555 if (atom->q != 0 || atom->qB != 0)
3563 if (EEL_FULL(ir->coulombtype))
3566 "You are using full electrostatics treatment %s for a system without charges.\n"
3567 "This costs a lot of performance for just processing zeros, consider using %s instead.\n",
3568 EELTYPE(ir->coulombtype), EELTYPE(eelCUT));
3569 warning(wi, err_buf);
3574 if (ir->coulombtype == eelCUT && ir->rcoulomb > 0 && !ir->implicit_solvent)
3577 "You are using a plain Coulomb cut-off, which might produce artifacts.\n"
3578 "You might want to consider using %s electrostatics.\n",
3580 warning_note(wi, err_buf);
3584 /* Generalized reaction field */
3585 if (ir->opts.ngtc == 0)
3587 sprintf(err_buf, "No temperature coupling while using coulombtype %s",
3589 CHECK(ir->coulombtype == eelGRF);
3593 sprintf(err_buf, "When using coulombtype = %s"
3594 " ref-t for temperature coupling should be > 0",
3596 CHECK((ir->coulombtype == eelGRF) && (ir->opts.ref_t[0] <= 0));
3599 if (ir->eI == eiSD1 &&
3600 (gmx_mtop_ftype_count(sys, F_CONSTR) > 0 ||
3601 gmx_mtop_ftype_count(sys, F_SETTLE) > 0))
3603 sprintf(warn_buf, "With constraints integrator %s is less accurate, consider using %s instead", ei_names[ir->eI], ei_names[eiSD2]);
3604 warning_note(wi, warn_buf);
3608 for (i = 0; (i < sys->groups.grps[egcACC].nr); i++)
3610 for (m = 0; (m < DIM); m++)
3612 if (fabs(ir->opts.acc[i][m]) > 1e-6)
3621 snew(mgrp, sys->groups.grps[egcACC].nr);
3622 aloop = gmx_mtop_atomloop_all_init(sys);
3623 while (gmx_mtop_atomloop_all_next(aloop, &i, &atom))
3625 mgrp[ggrpnr(&sys->groups, egcACC, i)] += atom->m;
3628 for (i = 0; (i < sys->groups.grps[egcACC].nr); i++)
3630 for (m = 0; (m < DIM); m++)
3632 acc[m] += ir->opts.acc[i][m]*mgrp[i];
3636 for (m = 0; (m < DIM); m++)
3638 if (fabs(acc[m]) > 1e-6)
3640 const char *dim[DIM] = { "X", "Y", "Z" };
3642 "Net Acceleration in %s direction, will %s be corrected\n",
3643 dim[m], ir->nstcomm != 0 ? "" : "not");
3644 if (ir->nstcomm != 0 && m < ndof_com(ir))
3647 for (i = 0; (i < sys->groups.grps[egcACC].nr); i++)
3649 ir->opts.acc[i][m] -= acc[m];
3657 if (ir->efep != efepNO && ir->fepvals->sc_alpha != 0 &&
3658 !gmx_within_tol(sys->ffparams.reppow, 12.0, 10*GMX_DOUBLE_EPS))
3660 gmx_fatal(FARGS, "Soft-core interactions are only supported with VdW repulsion power 12");
3663 if (ir->ePull != epullNO)
3665 if (ir->pull->grp[0].nat == 0)
3667 absolute_reference(ir, sys, FALSE, AbsRef);
3668 for (m = 0; m < DIM; m++)
3670 if (ir->pull->dim[m] && !AbsRef[m])
3672 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.");
3678 if (ir->pull->eGeom == epullgDIRPBC)
3680 for (i = 0; i < 3; i++)
3682 for (m = 0; m <= i; m++)
3684 if ((ir->epc != epcNO && ir->compress[i][m] != 0) ||
3685 ir->deform[i][m] != 0)
3687 for (g = 1; g < ir->pull->ngrp; g++)
3689 if (ir->pull->grp[g].vec[m] != 0)
3691 gmx_fatal(FARGS, "Can not have dynamic box while using pull geometry '%s' (dim %c)", EPULLGEOM(ir->pull->eGeom), 'x'+m);
3703 void double_check(t_inputrec *ir, matrix box, gmx_bool bConstr, warninp_t wi)
3707 char warn_buf[STRLEN];
3710 ptr = check_box(ir->ePBC, box);
3713 warning_error(wi, ptr);
3716 if (bConstr && ir->eConstrAlg == econtSHAKE)
3718 if (ir->shake_tol <= 0.0)
3720 sprintf(warn_buf, "ERROR: shake-tol must be > 0 instead of %g\n",
3722 warning_error(wi, warn_buf);
3725 if (IR_TWINRANGE(*ir) && ir->nstlist > 1)
3727 sprintf(warn_buf, "With twin-range cut-off's and SHAKE the virial and the pressure are incorrect.");
3728 if (ir->epc == epcNO)
3730 warning(wi, warn_buf);
3734 warning_error(wi, warn_buf);
3739 if ( (ir->eConstrAlg == econtLINCS) && bConstr)
3741 /* If we have Lincs constraints: */
3742 if (ir->eI == eiMD && ir->etc == etcNO &&
3743 ir->eConstrAlg == econtLINCS && ir->nLincsIter == 1)
3745 sprintf(warn_buf, "For energy conservation with LINCS, lincs_iter should be 2 or larger.\n");
3746 warning_note(wi, warn_buf);
3749 if ((ir->eI == eiCG || ir->eI == eiLBFGS) && (ir->nProjOrder < 8))
3751 sprintf(warn_buf, "For accurate %s with LINCS constraints, lincs-order should be 8 or more.", ei_names[ir->eI]);
3752 warning_note(wi, warn_buf);
3754 if (ir->epc == epcMTTK)
3756 warning_error(wi, "MTTK not compatible with lincs -- use shake instead.");
3760 if (ir->LincsWarnAngle > 90.0)
3762 sprintf(warn_buf, "lincs-warnangle can not be larger than 90 degrees, setting it to 90.\n");
3763 warning(wi, warn_buf);
3764 ir->LincsWarnAngle = 90.0;
3767 if (ir->ePBC != epbcNONE)
3769 if (ir->nstlist == 0)
3771 warning(wi, "With nstlist=0 atoms are only put into the box at step 0, therefore drifting atoms might cause the simulation to crash.");
3773 bTWIN = (ir->rlistlong > ir->rlist);
3774 if (ir->ns_type == ensGRID)
3776 if (sqr(ir->rlistlong) >= max_cutoff2(ir->ePBC, box))
3778 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",
3779 bTWIN ? (ir->rcoulomb == ir->rlistlong ? "rcoulomb" : "rvdw") : "rlist");
3780 warning_error(wi, warn_buf);
3785 min_size = min(box[XX][XX], min(box[YY][YY], box[ZZ][ZZ]));
3786 if (2*ir->rlistlong >= min_size)
3788 sprintf(warn_buf, "ERROR: One of the box lengths is smaller than twice the cut-off length. Increase the box size or decrease rlist.");
3789 warning_error(wi, warn_buf);
3792 fprintf(stderr, "Grid search might allow larger cut-off's than simple search with triclinic boxes.");
3799 void check_chargegroup_radii(const gmx_mtop_t *mtop, const t_inputrec *ir,
3803 real rvdw1, rvdw2, rcoul1, rcoul2;
3804 char warn_buf[STRLEN];
3806 calc_chargegroup_radii(mtop, x, &rvdw1, &rvdw2, &rcoul1, &rcoul2);
3810 printf("Largest charge group radii for Van der Waals: %5.3f, %5.3f nm\n",
3815 printf("Largest charge group radii for Coulomb: %5.3f, %5.3f nm\n",
3821 if (rvdw1 + rvdw2 > ir->rlist ||
3822 rcoul1 + rcoul2 > ir->rlist)
3824 sprintf(warn_buf, "The sum of the two largest charge group radii (%f) is larger than rlist (%f)\n", max(rvdw1+rvdw2, rcoul1+rcoul2), ir->rlist);
3825 warning(wi, warn_buf);
3829 /* Here we do not use the zero at cut-off macro,
3830 * since user defined interactions might purposely
3831 * not be zero at the cut-off.
3833 if (EVDW_IS_ZERO_AT_CUTOFF(ir->vdwtype) &&
3834 rvdw1 + rvdw2 > ir->rlist - ir->rvdw)
3836 sprintf(warn_buf, "The sum of the two largest charge group radii (%f) is larger than rlist (%f) - rvdw (%f)\n",
3838 ir->rlist, ir->rvdw);
3841 warning(wi, warn_buf);
3845 warning_note(wi, warn_buf);
3848 if (EEL_IS_ZERO_AT_CUTOFF(ir->coulombtype) &&
3849 rcoul1 + rcoul2 > ir->rlistlong - ir->rcoulomb)
3851 sprintf(warn_buf, "The sum of the two largest charge group radii (%f) is larger than %s (%f) - rcoulomb (%f)\n",
3853 ir->rlistlong > ir->rlist ? "rlistlong" : "rlist",
3854 ir->rlistlong, ir->rcoulomb);
3857 warning(wi, warn_buf);
3861 warning_note(wi, warn_buf);