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 (ir->coulombtype == eelPMESWITCH)
1081 if (ir->rcoulomb_switch/ir->rcoulomb < 0.9499)
1083 sprintf(warn_buf, "The switching range for %s should be 5%% or less, energy conservation will be good anyhow, since ewald_rtol = %g",
1084 eel_names[ir->coulombtype],
1086 warning(wi, warn_buf);
1090 if (EEL_FULL(ir->coulombtype))
1092 if (ir->coulombtype == eelPMESWITCH || ir->coulombtype == eelPMEUSER ||
1093 ir->coulombtype == eelPMEUSERSWITCH)
1095 sprintf(err_buf, "With coulombtype = %s, rcoulomb must be <= rlist",
1096 eel_names[ir->coulombtype]);
1097 CHECK(ir->rcoulomb > ir->rlist);
1099 else if (ir->cutoff_scheme == ecutsGROUP && ir->coulomb_modifier == eintmodNONE)
1101 if (ir->coulombtype == eelPME || ir->coulombtype == eelP3M_AD)
1104 "With coulombtype = %s (without modifier), rcoulomb must be equal to rlist,\n"
1105 "or rlistlong if nstcalclr=1. For optimal energy conservation,consider using\n"
1106 "a potential modifier.", eel_names[ir->coulombtype]);
1107 if (ir->nstcalclr == 1)
1109 CHECK(ir->rcoulomb != ir->rlist && ir->rcoulomb != ir->rlistlong);
1113 CHECK(ir->rcoulomb != ir->rlist);
1119 if (EEL_PME(ir->coulombtype))
1121 if (ir->pme_order < 3)
1123 warning_error(wi, "pme-order can not be smaller than 3");
1127 if (ir->nwall == 2 && EEL_FULL(ir->coulombtype))
1129 if (ir->ewald_geometry == eewg3D)
1131 sprintf(warn_buf, "With pbc=%s you should use ewald-geometry=%s",
1132 epbc_names[ir->ePBC], eewg_names[eewg3DC]);
1133 warning(wi, warn_buf);
1135 /* This check avoids extra pbc coding for exclusion corrections */
1136 sprintf(err_buf, "wall-ewald-zfac should be >= 2");
1137 CHECK(ir->wall_ewald_zfac < 2);
1140 if (EVDW_SWITCHED(ir->vdwtype))
1142 sprintf(err_buf, "With vdwtype = %s rvdw-switch must be < rvdw. Or, better - use a potential modifier.",
1143 evdw_names[ir->vdwtype]);
1144 CHECK(ir->rvdw_switch >= ir->rvdw);
1146 else if (ir->vdwtype == evdwCUT)
1148 if (ir->cutoff_scheme == ecutsGROUP && ir->vdw_modifier == eintmodNONE)
1150 sprintf(err_buf, "With vdwtype = %s, rvdw must be >= rlist unless you use a potential modifier", evdw_names[ir->vdwtype]);
1151 CHECK(ir->rlist > ir->rvdw);
1154 if (ir->cutoff_scheme == ecutsGROUP)
1156 if (((ir->coulomb_modifier != eintmodNONE && ir->rcoulomb == ir->rlist) ||
1157 (ir->vdw_modifier != eintmodNONE && ir->rvdw == ir->rlist)) &&
1160 warning_note(wi, "With exact cut-offs, rlist should be "
1161 "larger than rcoulomb and rvdw, so that there "
1162 "is a buffer region for particle motion "
1163 "between neighborsearch steps");
1166 if (EEL_IS_ZERO_AT_CUTOFF(ir->coulombtype)
1167 && (ir->rlistlong <= ir->rcoulomb))
1169 sprintf(warn_buf, "For energy conservation with switch/shift potentials, %s should be 0.1 to 0.3 nm larger than rcoulomb.",
1170 IR_TWINRANGE(*ir) ? "rlistlong" : "rlist");
1171 warning_note(wi, warn_buf);
1173 if (EVDW_SWITCHED(ir->vdwtype) && (ir->rlistlong <= ir->rvdw))
1175 sprintf(warn_buf, "For energy conservation with switch/shift potentials, %s should be 0.1 to 0.3 nm larger than rvdw.",
1176 IR_TWINRANGE(*ir) ? "rlistlong" : "rlist");
1177 warning_note(wi, warn_buf);
1181 if (ir->vdwtype == evdwUSER && ir->eDispCorr != edispcNO)
1183 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.");
1186 if (ir->nstlist == -1)
1188 sprintf(err_buf, "With nstlist=-1 rvdw and rcoulomb should be smaller than rlist to account for diffusion and possibly charge-group radii");
1189 CHECK(ir->rvdw >= ir->rlist || ir->rcoulomb >= ir->rlist);
1191 sprintf(err_buf, "nstlist can not be smaller than -1");
1192 CHECK(ir->nstlist < -1);
1194 if (ir->eI == eiLBFGS && (ir->coulombtype == eelCUT || ir->vdwtype == evdwCUT)
1197 warning(wi, "For efficient BFGS minimization, use switch/shift/pme instead of cut-off.");
1200 if (ir->eI == eiLBFGS && ir->nbfgscorr <= 0)
1202 warning(wi, "Using L-BFGS with nbfgscorr<=0 just gets you steepest descent.");
1205 /* ENERGY CONSERVATION */
1206 if (ir_NVE(ir) && ir->cutoff_scheme == ecutsGROUP)
1208 if (!EVDW_MIGHT_BE_ZERO_AT_CUTOFF(ir->vdwtype) && ir->rvdw > 0 && ir->vdw_modifier == eintmodNONE)
1210 sprintf(warn_buf, "You are using a cut-off for VdW interactions with NVE, for good energy conservation use vdwtype = %s (possibly with DispCorr)",
1211 evdw_names[evdwSHIFT]);
1212 warning_note(wi, warn_buf);
1214 if (!EEL_MIGHT_BE_ZERO_AT_CUTOFF(ir->coulombtype) && ir->rcoulomb > 0 && ir->coulomb_modifier == eintmodNONE)
1216 sprintf(warn_buf, "You are using a cut-off for electrostatics with NVE, for good energy conservation use coulombtype = %s or %s",
1217 eel_names[eelPMESWITCH], eel_names[eelRF_ZERO]);
1218 warning_note(wi, warn_buf);
1222 /* IMPLICIT SOLVENT */
1223 if (ir->coulombtype == eelGB_NOTUSED)
1225 ir->coulombtype = eelCUT;
1226 ir->implicit_solvent = eisGBSA;
1227 fprintf(stderr, "Note: Old option for generalized born electrostatics given:\n"
1228 "Changing coulombtype from \"generalized-born\" to \"cut-off\" and instead\n"
1229 "setting implicit-solvent value to \"GBSA\" in input section.\n");
1232 if (ir->sa_algorithm == esaSTILL)
1234 sprintf(err_buf, "Still SA algorithm not available yet, use %s or %s instead\n", esa_names[esaAPPROX], esa_names[esaNO]);
1235 CHECK(ir->sa_algorithm == esaSTILL);
1238 if (ir->implicit_solvent == eisGBSA)
1240 sprintf(err_buf, "With GBSA implicit solvent, rgbradii must be equal to rlist.");
1241 CHECK(ir->rgbradii != ir->rlist);
1243 if (ir->coulombtype != eelCUT)
1245 sprintf(err_buf, "With GBSA, coulombtype must be equal to %s\n", eel_names[eelCUT]);
1246 CHECK(ir->coulombtype != eelCUT);
1248 if (ir->vdwtype != evdwCUT)
1250 sprintf(err_buf, "With GBSA, vdw-type must be equal to %s\n", evdw_names[evdwCUT]);
1251 CHECK(ir->vdwtype != evdwCUT);
1253 if (ir->nstgbradii < 1)
1255 sprintf(warn_buf, "Using GBSA with nstgbradii<1, setting nstgbradii=1");
1256 warning_note(wi, warn_buf);
1259 if (ir->sa_algorithm == esaNO)
1261 sprintf(warn_buf, "No SA (non-polar) calculation requested together with GB. Are you sure this is what you want?\n");
1262 warning_note(wi, warn_buf);
1264 if (ir->sa_surface_tension < 0 && ir->sa_algorithm != esaNO)
1266 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");
1267 warning_note(wi, warn_buf);
1269 if (ir->gb_algorithm == egbSTILL)
1271 ir->sa_surface_tension = 0.0049 * CAL2JOULE * 100;
1275 ir->sa_surface_tension = 0.0054 * CAL2JOULE * 100;
1278 if (ir->sa_surface_tension == 0 && ir->sa_algorithm != esaNO)
1280 sprintf(err_buf, "Surface tension set to 0 while SA-calculation requested\n");
1281 CHECK(ir->sa_surface_tension == 0 && ir->sa_algorithm != esaNO);
1288 if (ir->cutoff_scheme != ecutsGROUP)
1290 warning_error(wi, "AdresS simulation supports only cutoff-scheme=group");
1294 warning_error(wi, "AdresS simulation supports only stochastic dynamics");
1296 if (ir->epc != epcNO)
1298 warning_error(wi, "AdresS simulation does not support pressure coupling");
1300 if (EEL_FULL(ir->coulombtype))
1302 warning_error(wi, "AdresS simulation does not support long-range electrostatics");
1307 /* count the number of text elemets separated by whitespace in a string.
1308 str = the input string
1309 maxptr = the maximum number of allowed elements
1310 ptr = the output array of pointers to the first character of each element
1311 returns: the number of elements. */
1312 int str_nelem(const char *str, int maxptr, char *ptr[])
1317 copy0 = strdup(str);
1320 while (*copy != '\0')
1324 gmx_fatal(FARGS, "Too many groups on line: '%s' (max is %d)",
1332 while ((*copy != '\0') && !isspace(*copy))
1351 /* interpret a number of doubles from a string and put them in an array,
1352 after allocating space for them.
1353 str = the input string
1354 n = the (pre-allocated) number of doubles read
1355 r = the output array of doubles. */
1356 static void parse_n_real(char *str, int *n, real **r)
1361 *n = str_nelem(str, MAXPTR, ptr);
1364 for (i = 0; i < *n; i++)
1366 (*r)[i] = strtod(ptr[i], NULL);
1370 static void do_fep_params(t_inputrec *ir, char fep_lambda[][STRLEN], char weights[STRLEN])
1373 int i, j, max_n_lambda, nweights, nfep[efptNR];
1374 t_lambda *fep = ir->fepvals;
1375 t_expanded *expand = ir->expandedvals;
1376 real **count_fep_lambdas;
1377 gmx_bool bOneLambda = TRUE;
1379 snew(count_fep_lambdas, efptNR);
1381 /* FEP input processing */
1382 /* first, identify the number of lambda values for each type.
1383 All that are nonzero must have the same number */
1385 for (i = 0; i < efptNR; i++)
1387 parse_n_real(fep_lambda[i], &(nfep[i]), &(count_fep_lambdas[i]));
1390 /* now, determine the number of components. All must be either zero, or equal. */
1393 for (i = 0; i < efptNR; i++)
1395 if (nfep[i] > max_n_lambda)
1397 max_n_lambda = nfep[i]; /* here's a nonzero one. All of them
1398 must have the same number if its not zero.*/
1403 for (i = 0; i < efptNR; i++)
1407 ir->fepvals->separate_dvdl[i] = FALSE;
1409 else if (nfep[i] == max_n_lambda)
1411 if (i != efptTEMPERATURE) /* we treat this differently -- not really a reason to compute the derivative with
1412 respect to the temperature currently */
1414 ir->fepvals->separate_dvdl[i] = TRUE;
1419 gmx_fatal(FARGS, "Number of lambdas (%d) for FEP type %s not equal to number of other types (%d)",
1420 nfep[i], efpt_names[i], max_n_lambda);
1423 /* we don't print out dhdl if the temperature is changing, since we can't correctly define dhdl in this case */
1424 ir->fepvals->separate_dvdl[efptTEMPERATURE] = FALSE;
1426 /* the number of lambdas is the number we've read in, which is either zero
1427 or the same for all */
1428 fep->n_lambda = max_n_lambda;
1430 /* allocate space for the array of lambda values */
1431 snew(fep->all_lambda, efptNR);
1432 /* if init_lambda is defined, we need to set lambda */
1433 if ((fep->init_lambda > 0) && (fep->n_lambda == 0))
1435 ir->fepvals->separate_dvdl[efptFEP] = TRUE;
1437 /* otherwise allocate the space for all of the lambdas, and transfer the data */
1438 for (i = 0; i < efptNR; i++)
1440 snew(fep->all_lambda[i], fep->n_lambda);
1441 if (nfep[i] > 0) /* if it's zero, then the count_fep_lambda arrays
1444 for (j = 0; j < fep->n_lambda; j++)
1446 fep->all_lambda[i][j] = (double)count_fep_lambdas[i][j];
1448 sfree(count_fep_lambdas[i]);
1451 sfree(count_fep_lambdas);
1453 /* "fep-vals" is either zero or the full number. If zero, we'll need to define fep-lambdas for internal
1454 bookkeeping -- for now, init_lambda */
1456 if ((nfep[efptFEP] == 0) && (fep->init_lambda >= 0))
1458 for (i = 0; i < fep->n_lambda; i++)
1460 fep->all_lambda[efptFEP][i] = fep->init_lambda;
1464 /* check to see if only a single component lambda is defined, and soft core is defined.
1465 In this case, turn on coulomb soft core */
1467 if (max_n_lambda == 0)
1473 for (i = 0; i < efptNR; i++)
1475 if ((nfep[i] != 0) && (i != efptFEP))
1481 if ((bOneLambda) && (fep->sc_alpha > 0))
1483 fep->bScCoul = TRUE;
1486 /* Fill in the others with the efptFEP if they are not explicitly
1487 specified (i.e. nfep[i] == 0). This means if fep is not defined,
1488 they are all zero. */
1490 for (i = 0; i < efptNR; i++)
1492 if ((nfep[i] == 0) && (i != efptFEP))
1494 for (j = 0; j < fep->n_lambda; j++)
1496 fep->all_lambda[i][j] = fep->all_lambda[efptFEP][j];
1502 /* make it easier if sc_r_power = 48 by increasing it to the 4th power, to be in the right scale. */
1503 if (fep->sc_r_power == 48)
1505 if (fep->sc_alpha > 0.1)
1507 gmx_fatal(FARGS, "sc_alpha (%f) for sc_r_power = 48 should usually be between 0.001 and 0.004", fep->sc_alpha);
1511 expand = ir->expandedvals;
1512 /* now read in the weights */
1513 parse_n_real(weights, &nweights, &(expand->init_lambda_weights));
1516 expand->bInit_weights = FALSE;
1517 snew(expand->init_lambda_weights, fep->n_lambda); /* initialize to zero */
1519 else if (nweights != fep->n_lambda)
1521 gmx_fatal(FARGS, "Number of weights (%d) is not equal to number of lambda values (%d)",
1522 nweights, fep->n_lambda);
1526 expand->bInit_weights = TRUE;
1528 if ((expand->nstexpanded < 0) && (ir->efep != efepNO))
1530 expand->nstexpanded = fep->nstdhdl;
1531 /* if you don't specify nstexpanded when doing expanded ensemble free energy calcs, it is set to nstdhdl */
1533 if ((expand->nstexpanded < 0) && ir->bSimTemp)
1535 expand->nstexpanded = 2*(int)(ir->opts.tau_t[0]/ir->delta_t);
1536 /* if you don't specify nstexpanded when doing expanded ensemble simulated tempering, it is set to
1537 2*tau_t just to be careful so it's not to frequent */
1542 static void do_simtemp_params(t_inputrec *ir)
1545 snew(ir->simtempvals->temperatures, ir->fepvals->n_lambda);
1546 GetSimTemps(ir->fepvals->n_lambda, ir->simtempvals, ir->fepvals->all_lambda[efptTEMPERATURE]);
1551 static void do_wall_params(t_inputrec *ir,
1552 char *wall_atomtype, char *wall_density,
1556 char *names[MAXPTR];
1559 opts->wall_atomtype[0] = NULL;
1560 opts->wall_atomtype[1] = NULL;
1562 ir->wall_atomtype[0] = -1;
1563 ir->wall_atomtype[1] = -1;
1564 ir->wall_density[0] = 0;
1565 ir->wall_density[1] = 0;
1569 nstr = str_nelem(wall_atomtype, MAXPTR, names);
1570 if (nstr != ir->nwall)
1572 gmx_fatal(FARGS, "Expected %d elements for wall_atomtype, found %d",
1575 for (i = 0; i < ir->nwall; i++)
1577 opts->wall_atomtype[i] = strdup(names[i]);
1580 if (ir->wall_type == ewt93 || ir->wall_type == ewt104)
1582 nstr = str_nelem(wall_density, MAXPTR, names);
1583 if (nstr != ir->nwall)
1585 gmx_fatal(FARGS, "Expected %d elements for wall-density, found %d", ir->nwall, nstr);
1587 for (i = 0; i < ir->nwall; i++)
1589 sscanf(names[i], "%lf", &dbl);
1592 gmx_fatal(FARGS, "wall-density[%d] = %f\n", i, dbl);
1594 ir->wall_density[i] = dbl;
1600 static void add_wall_energrps(gmx_groups_t *groups, int nwall, t_symtab *symtab)
1608 srenew(groups->grpname, groups->ngrpname+nwall);
1609 grps = &(groups->grps[egcENER]);
1610 srenew(grps->nm_ind, grps->nr+nwall);
1611 for (i = 0; i < nwall; i++)
1613 sprintf(str, "wall%d", i);
1614 groups->grpname[groups->ngrpname] = put_symtab(symtab, str);
1615 grps->nm_ind[grps->nr++] = groups->ngrpname++;
1620 void read_expandedparams(int *ninp_p, t_inpfile **inp_p,
1621 t_expanded *expand, warninp_t wi)
1623 int ninp, nerror = 0;
1629 /* read expanded ensemble parameters */
1630 CCTYPE ("expanded ensemble variables");
1631 ITYPE ("nstexpanded", expand->nstexpanded, -1);
1632 EETYPE("lmc-stats", expand->elamstats, elamstats_names);
1633 EETYPE("lmc-move", expand->elmcmove, elmcmove_names);
1634 EETYPE("lmc-weights-equil", expand->elmceq, elmceq_names);
1635 ITYPE ("weight-equil-number-all-lambda", expand->equil_n_at_lam, -1);
1636 ITYPE ("weight-equil-number-samples", expand->equil_samples, -1);
1637 ITYPE ("weight-equil-number-steps", expand->equil_steps, -1);
1638 RTYPE ("weight-equil-wl-delta", expand->equil_wl_delta, -1);
1639 RTYPE ("weight-equil-count-ratio", expand->equil_ratio, -1);
1640 CCTYPE("Seed for Monte Carlo in lambda space");
1641 ITYPE ("lmc-seed", expand->lmc_seed, -1);
1642 RTYPE ("mc-temperature", expand->mc_temp, -1);
1643 ITYPE ("lmc-repeats", expand->lmc_repeats, 1);
1644 ITYPE ("lmc-gibbsdelta", expand->gibbsdeltalam, -1);
1645 ITYPE ("lmc-forced-nstart", expand->lmc_forced_nstart, 0);
1646 EETYPE("symmetrized-transition-matrix", expand->bSymmetrizedTMatrix, yesno_names);
1647 ITYPE("nst-transition-matrix", expand->nstTij, -1);
1648 ITYPE ("mininum-var-min", expand->minvarmin, 100); /*default is reasonable */
1649 ITYPE ("weight-c-range", expand->c_range, 0); /* default is just C=0 */
1650 RTYPE ("wl-scale", expand->wl_scale, 0.8);
1651 RTYPE ("wl-ratio", expand->wl_ratio, 0.8);
1652 RTYPE ("init-wl-delta", expand->init_wl_delta, 1.0);
1653 EETYPE("wl-oneovert", expand->bWLoneovert, yesno_names);
1661 void get_ir(const char *mdparin, const char *mdparout,
1662 t_inputrec *ir, t_gromppopts *opts,
1666 double dumdub[2][6];
1670 char warn_buf[STRLEN];
1671 t_lambda *fep = ir->fepvals;
1672 t_expanded *expand = ir->expandedvals;
1674 inp = read_inpfile(mdparin, &ninp, NULL, wi);
1676 snew(dumstr[0], STRLEN);
1677 snew(dumstr[1], STRLEN);
1679 /* remove the following deprecated commands */
1682 REM_TYPE("domain-decomposition");
1683 REM_TYPE("andersen-seed");
1685 REM_TYPE("dihre-fc");
1686 REM_TYPE("dihre-tau");
1687 REM_TYPE("nstdihreout");
1688 REM_TYPE("nstcheckpoint");
1690 /* replace the following commands with the clearer new versions*/
1691 REPL_TYPE("unconstrained-start", "continuation");
1692 REPL_TYPE("foreign-lambda", "fep-lambdas");
1694 CCTYPE ("VARIOUS PREPROCESSING OPTIONS");
1695 CTYPE ("Preprocessor information: use cpp syntax.");
1696 CTYPE ("e.g.: -I/home/joe/doe -I/home/mary/roe");
1697 STYPE ("include", opts->include, NULL);
1698 CTYPE ("e.g.: -DPOSRES -DFLEXIBLE (note these variable names are case sensitive)");
1699 STYPE ("define", opts->define, NULL);
1701 CCTYPE ("RUN CONTROL PARAMETERS");
1702 EETYPE("integrator", ir->eI, ei_names);
1703 CTYPE ("Start time and timestep in ps");
1704 RTYPE ("tinit", ir->init_t, 0.0);
1705 RTYPE ("dt", ir->delta_t, 0.001);
1706 STEPTYPE ("nsteps", ir->nsteps, 0);
1707 CTYPE ("For exact run continuation or redoing part of a run");
1708 STEPTYPE ("init-step", ir->init_step, 0);
1709 CTYPE ("Part index is updated automatically on checkpointing (keeps files separate)");
1710 ITYPE ("simulation-part", ir->simulation_part, 1);
1711 CTYPE ("mode for center of mass motion removal");
1712 EETYPE("comm-mode", ir->comm_mode, ecm_names);
1713 CTYPE ("number of steps for center of mass motion removal");
1714 ITYPE ("nstcomm", ir->nstcomm, 100);
1715 CTYPE ("group(s) for center of mass motion removal");
1716 STYPE ("comm-grps", vcm, NULL);
1718 CCTYPE ("LANGEVIN DYNAMICS OPTIONS");
1719 CTYPE ("Friction coefficient (amu/ps) and random seed");
1720 RTYPE ("bd-fric", ir->bd_fric, 0.0);
1721 ITYPE ("ld-seed", ir->ld_seed, 1993);
1724 CCTYPE ("ENERGY MINIMIZATION OPTIONS");
1725 CTYPE ("Force tolerance and initial step-size");
1726 RTYPE ("emtol", ir->em_tol, 10.0);
1727 RTYPE ("emstep", ir->em_stepsize, 0.01);
1728 CTYPE ("Max number of iterations in relax-shells");
1729 ITYPE ("niter", ir->niter, 20);
1730 CTYPE ("Step size (ps^2) for minimization of flexible constraints");
1731 RTYPE ("fcstep", ir->fc_stepsize, 0);
1732 CTYPE ("Frequency of steepest descents steps when doing CG");
1733 ITYPE ("nstcgsteep", ir->nstcgsteep, 1000);
1734 ITYPE ("nbfgscorr", ir->nbfgscorr, 10);
1736 CCTYPE ("TEST PARTICLE INSERTION OPTIONS");
1737 RTYPE ("rtpi", ir->rtpi, 0.05);
1739 /* Output options */
1740 CCTYPE ("OUTPUT CONTROL OPTIONS");
1741 CTYPE ("Output frequency for coords (x), velocities (v) and forces (f)");
1742 ITYPE ("nstxout", ir->nstxout, 0);
1743 ITYPE ("nstvout", ir->nstvout, 0);
1744 ITYPE ("nstfout", ir->nstfout, 0);
1745 ir->nstcheckpoint = 1000;
1746 CTYPE ("Output frequency for energies to log file and energy file");
1747 ITYPE ("nstlog", ir->nstlog, 1000);
1748 ITYPE ("nstcalcenergy", ir->nstcalcenergy, 100);
1749 ITYPE ("nstenergy", ir->nstenergy, 1000);
1750 CTYPE ("Output frequency and precision for .xtc file");
1751 ITYPE ("nstxtcout", ir->nstxtcout, 0);
1752 RTYPE ("xtc-precision", ir->xtcprec, 1000.0);
1753 CTYPE ("This selects the subset of atoms for the .xtc file. You can");
1754 CTYPE ("select multiple groups. By default all atoms will be written.");
1755 STYPE ("xtc-grps", xtc_grps, NULL);
1756 CTYPE ("Selection of energy groups");
1757 STYPE ("energygrps", energy, NULL);
1759 /* Neighbor searching */
1760 CCTYPE ("NEIGHBORSEARCHING PARAMETERS");
1761 CTYPE ("cut-off scheme (group: using charge groups, Verlet: particle based cut-offs)");
1762 EETYPE("cutoff-scheme", ir->cutoff_scheme, ecutscheme_names);
1763 CTYPE ("nblist update frequency");
1764 ITYPE ("nstlist", ir->nstlist, 10);
1765 CTYPE ("ns algorithm (simple or grid)");
1766 EETYPE("ns-type", ir->ns_type, ens_names);
1767 /* set ndelta to the optimal value of 2 */
1769 CTYPE ("Periodic boundary conditions: xyz, no, xy");
1770 EETYPE("pbc", ir->ePBC, epbc_names);
1771 EETYPE("periodic-molecules", ir->bPeriodicMols, yesno_names);
1772 CTYPE ("Allowed energy drift due to the Verlet buffer in kJ/mol/ps per atom,");
1773 CTYPE ("a value of -1 means: use rlist");
1774 RTYPE("verlet-buffer-drift", ir->verletbuf_drift, 0.005);
1775 CTYPE ("nblist cut-off");
1776 RTYPE ("rlist", ir->rlist, 1.0);
1777 CTYPE ("long-range cut-off for switched potentials");
1778 RTYPE ("rlistlong", ir->rlistlong, -1);
1779 ITYPE ("nstcalclr", ir->nstcalclr, -1);
1781 /* Electrostatics */
1782 CCTYPE ("OPTIONS FOR ELECTROSTATICS AND VDW");
1783 CTYPE ("Method for doing electrostatics");
1784 EETYPE("coulombtype", ir->coulombtype, eel_names);
1785 EETYPE("coulomb-modifier", ir->coulomb_modifier, eintmod_names);
1786 CTYPE ("cut-off lengths");
1787 RTYPE ("rcoulomb-switch", ir->rcoulomb_switch, 0.0);
1788 RTYPE ("rcoulomb", ir->rcoulomb, 1.0);
1789 CTYPE ("Relative dielectric constant for the medium and the reaction field");
1790 RTYPE ("epsilon-r", ir->epsilon_r, 1.0);
1791 RTYPE ("epsilon-rf", ir->epsilon_rf, 0.0);
1792 CTYPE ("Method for doing Van der Waals");
1793 EETYPE("vdw-type", ir->vdwtype, evdw_names);
1794 EETYPE("vdw-modifier", ir->vdw_modifier, eintmod_names);
1795 CTYPE ("cut-off lengths");
1796 RTYPE ("rvdw-switch", ir->rvdw_switch, 0.0);
1797 RTYPE ("rvdw", ir->rvdw, 1.0);
1798 CTYPE ("Apply long range dispersion corrections for Energy and Pressure");
1799 EETYPE("DispCorr", ir->eDispCorr, edispc_names);
1800 CTYPE ("Extension of the potential lookup tables beyond the cut-off");
1801 RTYPE ("table-extension", ir->tabext, 1.0);
1802 CTYPE ("Separate tables between energy group pairs");
1803 STYPE ("energygrp-table", egptable, NULL);
1804 CTYPE ("Spacing for the PME/PPPM FFT grid");
1805 RTYPE ("fourierspacing", ir->fourier_spacing, 0.12);
1806 CTYPE ("FFT grid size, when a value is 0 fourierspacing will be used");
1807 ITYPE ("fourier-nx", ir->nkx, 0);
1808 ITYPE ("fourier-ny", ir->nky, 0);
1809 ITYPE ("fourier-nz", ir->nkz, 0);
1810 CTYPE ("EWALD/PME/PPPM parameters");
1811 ITYPE ("pme-order", ir->pme_order, 4);
1812 RTYPE ("ewald-rtol", ir->ewald_rtol, 0.00001);
1813 EETYPE("ewald-geometry", ir->ewald_geometry, eewg_names);
1814 RTYPE ("epsilon-surface", ir->epsilon_surface, 0.0);
1815 EETYPE("optimize-fft", ir->bOptFFT, yesno_names);
1817 CCTYPE("IMPLICIT SOLVENT ALGORITHM");
1818 EETYPE("implicit-solvent", ir->implicit_solvent, eis_names);
1820 CCTYPE ("GENERALIZED BORN ELECTROSTATICS");
1821 CTYPE ("Algorithm for calculating Born radii");
1822 EETYPE("gb-algorithm", ir->gb_algorithm, egb_names);
1823 CTYPE ("Frequency of calculating the Born radii inside rlist");
1824 ITYPE ("nstgbradii", ir->nstgbradii, 1);
1825 CTYPE ("Cutoff for Born radii calculation; the contribution from atoms");
1826 CTYPE ("between rlist and rgbradii is updated every nstlist steps");
1827 RTYPE ("rgbradii", ir->rgbradii, 1.0);
1828 CTYPE ("Dielectric coefficient of the implicit solvent");
1829 RTYPE ("gb-epsilon-solvent", ir->gb_epsilon_solvent, 80.0);
1830 CTYPE ("Salt concentration in M for Generalized Born models");
1831 RTYPE ("gb-saltconc", ir->gb_saltconc, 0.0);
1832 CTYPE ("Scaling factors used in the OBC GB model. Default values are OBC(II)");
1833 RTYPE ("gb-obc-alpha", ir->gb_obc_alpha, 1.0);
1834 RTYPE ("gb-obc-beta", ir->gb_obc_beta, 0.8);
1835 RTYPE ("gb-obc-gamma", ir->gb_obc_gamma, 4.85);
1836 RTYPE ("gb-dielectric-offset", ir->gb_dielectric_offset, 0.009);
1837 EETYPE("sa-algorithm", ir->sa_algorithm, esa_names);
1838 CTYPE ("Surface tension (kJ/mol/nm^2) for the SA (nonpolar surface) part of GBSA");
1839 CTYPE ("The value -1 will set default value for Still/HCT/OBC GB-models.");
1840 RTYPE ("sa-surface-tension", ir->sa_surface_tension, -1);
1842 /* Coupling stuff */
1843 CCTYPE ("OPTIONS FOR WEAK COUPLING ALGORITHMS");
1844 CTYPE ("Temperature coupling");
1845 EETYPE("tcoupl", ir->etc, etcoupl_names);
1846 ITYPE ("nsttcouple", ir->nsttcouple, -1);
1847 ITYPE("nh-chain-length", ir->opts.nhchainlength, NHCHAINLENGTH);
1848 EETYPE("print-nose-hoover-chain-variables", ir->bPrintNHChains, yesno_names);
1849 CTYPE ("Groups to couple separately");
1850 STYPE ("tc-grps", tcgrps, NULL);
1851 CTYPE ("Time constant (ps) and reference temperature (K)");
1852 STYPE ("tau-t", tau_t, NULL);
1853 STYPE ("ref-t", ref_t, NULL);
1854 CTYPE ("pressure coupling");
1855 EETYPE("pcoupl", ir->epc, epcoupl_names);
1856 EETYPE("pcoupltype", ir->epct, epcoupltype_names);
1857 ITYPE ("nstpcouple", ir->nstpcouple, -1);
1858 CTYPE ("Time constant (ps), compressibility (1/bar) and reference P (bar)");
1859 RTYPE ("tau-p", ir->tau_p, 1.0);
1860 STYPE ("compressibility", dumstr[0], NULL);
1861 STYPE ("ref-p", dumstr[1], NULL);
1862 CTYPE ("Scaling of reference coordinates, No, All or COM");
1863 EETYPE ("refcoord-scaling", ir->refcoord_scaling, erefscaling_names);
1866 CCTYPE ("OPTIONS FOR QMMM calculations");
1867 EETYPE("QMMM", ir->bQMMM, yesno_names);
1868 CTYPE ("Groups treated Quantum Mechanically");
1869 STYPE ("QMMM-grps", QMMM, NULL);
1870 CTYPE ("QM method");
1871 STYPE("QMmethod", QMmethod, NULL);
1872 CTYPE ("QMMM scheme");
1873 EETYPE("QMMMscheme", ir->QMMMscheme, eQMMMscheme_names);
1874 CTYPE ("QM basisset");
1875 STYPE("QMbasis", QMbasis, NULL);
1876 CTYPE ("QM charge");
1877 STYPE ("QMcharge", QMcharge, NULL);
1878 CTYPE ("QM multiplicity");
1879 STYPE ("QMmult", QMmult, NULL);
1880 CTYPE ("Surface Hopping");
1881 STYPE ("SH", bSH, NULL);
1882 CTYPE ("CAS space options");
1883 STYPE ("CASorbitals", CASorbitals, NULL);
1884 STYPE ("CASelectrons", CASelectrons, NULL);
1885 STYPE ("SAon", SAon, NULL);
1886 STYPE ("SAoff", SAoff, NULL);
1887 STYPE ("SAsteps", SAsteps, NULL);
1888 CTYPE ("Scale factor for MM charges");
1889 RTYPE ("MMChargeScaleFactor", ir->scalefactor, 1.0);
1890 CTYPE ("Optimization of QM subsystem");
1891 STYPE ("bOPT", bOPT, NULL);
1892 STYPE ("bTS", bTS, NULL);
1894 /* Simulated annealing */
1895 CCTYPE("SIMULATED ANNEALING");
1896 CTYPE ("Type of annealing for each temperature group (no/single/periodic)");
1897 STYPE ("annealing", anneal, NULL);
1898 CTYPE ("Number of time points to use for specifying annealing in each group");
1899 STYPE ("annealing-npoints", anneal_npoints, NULL);
1900 CTYPE ("List of times at the annealing points for each group");
1901 STYPE ("annealing-time", anneal_time, NULL);
1902 CTYPE ("Temp. at each annealing point, for each group.");
1903 STYPE ("annealing-temp", anneal_temp, NULL);
1906 CCTYPE ("GENERATE VELOCITIES FOR STARTUP RUN");
1907 EETYPE("gen-vel", opts->bGenVel, yesno_names);
1908 RTYPE ("gen-temp", opts->tempi, 300.0);
1909 ITYPE ("gen-seed", opts->seed, 173529);
1912 CCTYPE ("OPTIONS FOR BONDS");
1913 EETYPE("constraints", opts->nshake, constraints);
1914 CTYPE ("Type of constraint algorithm");
1915 EETYPE("constraint-algorithm", ir->eConstrAlg, econstr_names);
1916 CTYPE ("Do not constrain the start configuration");
1917 EETYPE("continuation", ir->bContinuation, yesno_names);
1918 CTYPE ("Use successive overrelaxation to reduce the number of shake iterations");
1919 EETYPE("Shake-SOR", ir->bShakeSOR, yesno_names);
1920 CTYPE ("Relative tolerance of shake");
1921 RTYPE ("shake-tol", ir->shake_tol, 0.0001);
1922 CTYPE ("Highest order in the expansion of the constraint coupling matrix");
1923 ITYPE ("lincs-order", ir->nProjOrder, 4);
1924 CTYPE ("Number of iterations in the final step of LINCS. 1 is fine for");
1925 CTYPE ("normal simulations, but use 2 to conserve energy in NVE runs.");
1926 CTYPE ("For energy minimization with constraints it should be 4 to 8.");
1927 ITYPE ("lincs-iter", ir->nLincsIter, 1);
1928 CTYPE ("Lincs will write a warning to the stderr if in one step a bond");
1929 CTYPE ("rotates over more degrees than");
1930 RTYPE ("lincs-warnangle", ir->LincsWarnAngle, 30.0);
1931 CTYPE ("Convert harmonic bonds to morse potentials");
1932 EETYPE("morse", opts->bMorse, yesno_names);
1934 /* Energy group exclusions */
1935 CCTYPE ("ENERGY GROUP EXCLUSIONS");
1936 CTYPE ("Pairs of energy groups for which all non-bonded interactions are excluded");
1937 STYPE ("energygrp-excl", egpexcl, NULL);
1941 CTYPE ("Number of walls, type, atom types, densities and box-z scale factor for Ewald");
1942 ITYPE ("nwall", ir->nwall, 0);
1943 EETYPE("wall-type", ir->wall_type, ewt_names);
1944 RTYPE ("wall-r-linpot", ir->wall_r_linpot, -1);
1945 STYPE ("wall-atomtype", wall_atomtype, NULL);
1946 STYPE ("wall-density", wall_density, NULL);
1947 RTYPE ("wall-ewald-zfac", ir->wall_ewald_zfac, 3);
1950 CCTYPE("COM PULLING");
1951 CTYPE("Pull type: no, umbrella, constraint or constant-force");
1952 EETYPE("pull", ir->ePull, epull_names);
1953 if (ir->ePull != epullNO)
1956 pull_grp = read_pullparams(&ninp, &inp, ir->pull, &opts->pull_start, wi);
1959 /* Enforced rotation */
1960 CCTYPE("ENFORCED ROTATION");
1961 CTYPE("Enforced rotation: No or Yes");
1962 EETYPE("rotation", ir->bRot, yesno_names);
1966 rot_grp = read_rotparams(&ninp, &inp, ir->rot, wi);
1970 CCTYPE("NMR refinement stuff");
1971 CTYPE ("Distance restraints type: No, Simple or Ensemble");
1972 EETYPE("disre", ir->eDisre, edisre_names);
1973 CTYPE ("Force weighting of pairs in one distance restraint: Conservative or Equal");
1974 EETYPE("disre-weighting", ir->eDisreWeighting, edisreweighting_names);
1975 CTYPE ("Use sqrt of the time averaged times the instantaneous violation");
1976 EETYPE("disre-mixed", ir->bDisreMixed, yesno_names);
1977 RTYPE ("disre-fc", ir->dr_fc, 1000.0);
1978 RTYPE ("disre-tau", ir->dr_tau, 0.0);
1979 CTYPE ("Output frequency for pair distances to energy file");
1980 ITYPE ("nstdisreout", ir->nstdisreout, 100);
1981 CTYPE ("Orientation restraints: No or Yes");
1982 EETYPE("orire", opts->bOrire, yesno_names);
1983 CTYPE ("Orientation restraints force constant and tau for time averaging");
1984 RTYPE ("orire-fc", ir->orires_fc, 0.0);
1985 RTYPE ("orire-tau", ir->orires_tau, 0.0);
1986 STYPE ("orire-fitgrp", orirefitgrp, NULL);
1987 CTYPE ("Output frequency for trace(SD) and S to energy file");
1988 ITYPE ("nstorireout", ir->nstorireout, 100);
1990 /* free energy variables */
1991 CCTYPE ("Free energy variables");
1992 EETYPE("free-energy", ir->efep, efep_names);
1993 STYPE ("couple-moltype", couple_moltype, NULL);
1994 EETYPE("couple-lambda0", opts->couple_lam0, couple_lam);
1995 EETYPE("couple-lambda1", opts->couple_lam1, couple_lam);
1996 EETYPE("couple-intramol", opts->bCoupleIntra, yesno_names);
1998 RTYPE ("init-lambda", fep->init_lambda, -1); /* start with -1 so
2000 it was not entered */
2001 ITYPE ("init-lambda-state", fep->init_fep_state, -1);
2002 RTYPE ("delta-lambda", fep->delta_lambda, 0.0);
2003 ITYPE ("nstdhdl", fep->nstdhdl, 50);
2004 STYPE ("fep-lambdas", fep_lambda[efptFEP], NULL);
2005 STYPE ("mass-lambdas", fep_lambda[efptMASS], NULL);
2006 STYPE ("coul-lambdas", fep_lambda[efptCOUL], NULL);
2007 STYPE ("vdw-lambdas", fep_lambda[efptVDW], NULL);
2008 STYPE ("bonded-lambdas", fep_lambda[efptBONDED], NULL);
2009 STYPE ("restraint-lambdas", fep_lambda[efptRESTRAINT], NULL);
2010 STYPE ("temperature-lambdas", fep_lambda[efptTEMPERATURE], NULL);
2011 ITYPE ("calc-lambda-neighbors", fep->lambda_neighbors, 1);
2012 STYPE ("init-lambda-weights", lambda_weights, NULL);
2013 EETYPE("dhdl-print-energy", fep->bPrintEnergy, yesno_names);
2014 RTYPE ("sc-alpha", fep->sc_alpha, 0.0);
2015 ITYPE ("sc-power", fep->sc_power, 1);
2016 RTYPE ("sc-r-power", fep->sc_r_power, 6.0);
2017 RTYPE ("sc-sigma", fep->sc_sigma, 0.3);
2018 EETYPE("sc-coul", fep->bScCoul, yesno_names);
2019 ITYPE ("dh_hist_size", fep->dh_hist_size, 0);
2020 RTYPE ("dh_hist_spacing", fep->dh_hist_spacing, 0.1);
2021 EETYPE("separate-dhdl-file", fep->separate_dhdl_file,
2022 separate_dhdl_file_names);
2023 EETYPE("dhdl-derivatives", fep->dhdl_derivatives, dhdl_derivatives_names);
2024 ITYPE ("dh_hist_size", fep->dh_hist_size, 0);
2025 RTYPE ("dh_hist_spacing", fep->dh_hist_spacing, 0.1);
2027 /* Non-equilibrium MD stuff */
2028 CCTYPE("Non-equilibrium MD stuff");
2029 STYPE ("acc-grps", accgrps, NULL);
2030 STYPE ("accelerate", acc, NULL);
2031 STYPE ("freezegrps", freeze, NULL);
2032 STYPE ("freezedim", frdim, NULL);
2033 RTYPE ("cos-acceleration", ir->cos_accel, 0);
2034 STYPE ("deform", deform, NULL);
2036 /* simulated tempering variables */
2037 CCTYPE("simulated tempering variables");
2038 EETYPE("simulated-tempering", ir->bSimTemp, yesno_names);
2039 EETYPE("simulated-tempering-scaling", ir->simtempvals->eSimTempScale, esimtemp_names);
2040 RTYPE("sim-temp-low", ir->simtempvals->simtemp_low, 300.0);
2041 RTYPE("sim-temp-high", ir->simtempvals->simtemp_high, 300.0);
2043 /* expanded ensemble variables */
2044 if (ir->efep == efepEXPANDED || ir->bSimTemp)
2046 read_expandedparams(&ninp, &inp, expand, wi);
2049 /* Electric fields */
2050 CCTYPE("Electric fields");
2051 CTYPE ("Format is number of terms (int) and for all terms an amplitude (real)");
2052 CTYPE ("and a phase angle (real)");
2053 STYPE ("E-x", efield_x, NULL);
2054 STYPE ("E-xt", efield_xt, NULL);
2055 STYPE ("E-y", efield_y, NULL);
2056 STYPE ("E-yt", efield_yt, NULL);
2057 STYPE ("E-z", efield_z, NULL);
2058 STYPE ("E-zt", efield_zt, NULL);
2060 /* AdResS defined thingies */
2061 CCTYPE ("AdResS parameters");
2062 EETYPE("adress", ir->bAdress, yesno_names);
2065 snew(ir->adress, 1);
2066 read_adressparams(&ninp, &inp, ir->adress, wi);
2069 /* User defined thingies */
2070 CCTYPE ("User defined thingies");
2071 STYPE ("user1-grps", user1, NULL);
2072 STYPE ("user2-grps", user2, NULL);
2073 ITYPE ("userint1", ir->userint1, 0);
2074 ITYPE ("userint2", ir->userint2, 0);
2075 ITYPE ("userint3", ir->userint3, 0);
2076 ITYPE ("userint4", ir->userint4, 0);
2077 RTYPE ("userreal1", ir->userreal1, 0);
2078 RTYPE ("userreal2", ir->userreal2, 0);
2079 RTYPE ("userreal3", ir->userreal3, 0);
2080 RTYPE ("userreal4", ir->userreal4, 0);
2083 write_inpfile(mdparout, ninp, inp, FALSE, wi);
2084 for (i = 0; (i < ninp); i++)
2087 sfree(inp[i].value);
2091 /* Process options if necessary */
2092 for (m = 0; m < 2; m++)
2094 for (i = 0; i < 2*DIM; i++)
2103 if (sscanf(dumstr[m], "%lf", &(dumdub[m][XX])) != 1)
2105 warning_error(wi, "Pressure coupling not enough values (I need 1)");
2107 dumdub[m][YY] = dumdub[m][ZZ] = dumdub[m][XX];
2109 case epctSEMIISOTROPIC:
2110 case epctSURFACETENSION:
2111 if (sscanf(dumstr[m], "%lf%lf",
2112 &(dumdub[m][XX]), &(dumdub[m][ZZ])) != 2)
2114 warning_error(wi, "Pressure coupling not enough values (I need 2)");
2116 dumdub[m][YY] = dumdub[m][XX];
2118 case epctANISOTROPIC:
2119 if (sscanf(dumstr[m], "%lf%lf%lf%lf%lf%lf",
2120 &(dumdub[m][XX]), &(dumdub[m][YY]), &(dumdub[m][ZZ]),
2121 &(dumdub[m][3]), &(dumdub[m][4]), &(dumdub[m][5])) != 6)
2123 warning_error(wi, "Pressure coupling not enough values (I need 6)");
2127 gmx_fatal(FARGS, "Pressure coupling type %s not implemented yet",
2128 epcoupltype_names[ir->epct]);
2132 clear_mat(ir->ref_p);
2133 clear_mat(ir->compress);
2134 for (i = 0; i < DIM; i++)
2136 ir->ref_p[i][i] = dumdub[1][i];
2137 ir->compress[i][i] = dumdub[0][i];
2139 if (ir->epct == epctANISOTROPIC)
2141 ir->ref_p[XX][YY] = dumdub[1][3];
2142 ir->ref_p[XX][ZZ] = dumdub[1][4];
2143 ir->ref_p[YY][ZZ] = dumdub[1][5];
2144 if (ir->ref_p[XX][YY] != 0 && ir->ref_p[XX][ZZ] != 0 && ir->ref_p[YY][ZZ] != 0)
2146 warning(wi, "All off-diagonal reference pressures are non-zero. Are you sure you want to apply a threefold shear stress?\n");
2148 ir->compress[XX][YY] = dumdub[0][3];
2149 ir->compress[XX][ZZ] = dumdub[0][4];
2150 ir->compress[YY][ZZ] = dumdub[0][5];
2151 for (i = 0; i < DIM; i++)
2153 for (m = 0; m < i; m++)
2155 ir->ref_p[i][m] = ir->ref_p[m][i];
2156 ir->compress[i][m] = ir->compress[m][i];
2161 if (ir->comm_mode == ecmNO)
2166 opts->couple_moltype = NULL;
2167 if (strlen(couple_moltype) > 0)
2169 if (ir->efep != efepNO)
2171 opts->couple_moltype = strdup(couple_moltype);
2172 if (opts->couple_lam0 == opts->couple_lam1)
2174 warning(wi, "The lambda=0 and lambda=1 states for coupling are identical");
2176 if (ir->eI == eiMD && (opts->couple_lam0 == ecouplamNONE ||
2177 opts->couple_lam1 == ecouplamNONE))
2179 warning(wi, "For proper sampling of the (nearly) decoupled state, stochastic dynamics should be used");
2184 warning(wi, "Can not couple a molecule with free_energy = no");
2187 /* FREE ENERGY AND EXPANDED ENSEMBLE OPTIONS */
2188 if (ir->efep != efepNO)
2190 if (fep->delta_lambda > 0)
2192 ir->efep = efepSLOWGROWTH;
2198 fep->bPrintEnergy = TRUE;
2199 /* always print out the energy to dhdl if we are doing expanded ensemble, since we need the total energy
2200 if the temperature is changing. */
2203 if ((ir->efep != efepNO) || ir->bSimTemp)
2205 ir->bExpanded = FALSE;
2206 if ((ir->efep == efepEXPANDED) || ir->bSimTemp)
2208 ir->bExpanded = TRUE;
2210 do_fep_params(ir, fep_lambda, lambda_weights);
2211 if (ir->bSimTemp) /* done after fep params */
2213 do_simtemp_params(ir);
2218 ir->fepvals->n_lambda = 0;
2221 /* WALL PARAMETERS */
2223 do_wall_params(ir, wall_atomtype, wall_density, opts);
2225 /* ORIENTATION RESTRAINT PARAMETERS */
2227 if (opts->bOrire && str_nelem(orirefitgrp, MAXPTR, NULL) != 1)
2229 warning_error(wi, "ERROR: Need one orientation restraint fit group\n");
2232 /* DEFORMATION PARAMETERS */
2234 clear_mat(ir->deform);
2235 for (i = 0; i < 6; i++)
2239 m = sscanf(deform, "%lf %lf %lf %lf %lf %lf",
2240 &(dumdub[0][0]), &(dumdub[0][1]), &(dumdub[0][2]),
2241 &(dumdub[0][3]), &(dumdub[0][4]), &(dumdub[0][5]));
2242 for (i = 0; i < 3; i++)
2244 ir->deform[i][i] = dumdub[0][i];
2246 ir->deform[YY][XX] = dumdub[0][3];
2247 ir->deform[ZZ][XX] = dumdub[0][4];
2248 ir->deform[ZZ][YY] = dumdub[0][5];
2249 if (ir->epc != epcNO)
2251 for (i = 0; i < 3; i++)
2253 for (j = 0; j <= i; j++)
2255 if (ir->deform[i][j] != 0 && ir->compress[i][j] != 0)
2257 warning_error(wi, "A box element has deform set and compressibility > 0");
2261 for (i = 0; i < 3; i++)
2263 for (j = 0; j < i; j++)
2265 if (ir->deform[i][j] != 0)
2267 for (m = j; m < DIM; m++)
2269 if (ir->compress[m][j] != 0)
2271 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.");
2272 warning(wi, warn_buf);
2284 static int search_QMstring(char *s, int ng, const char *gn[])
2286 /* same as normal search_string, but this one searches QM strings */
2289 for (i = 0; (i < ng); i++)
2291 if (gmx_strcasecmp(s, gn[i]) == 0)
2297 gmx_fatal(FARGS, "this QM method or basisset (%s) is not implemented\n!", s);
2301 } /* search_QMstring */
2304 int search_string(char *s, int ng, char *gn[])
2308 for (i = 0; (i < ng); i++)
2310 if (gmx_strcasecmp(s, gn[i]) == 0)
2317 "Group %s referenced in the .mdp file was not found in the index file.\n"
2318 "Group names must match either [moleculetype] names or custom index group\n"
2319 "names, in which case you must supply an index file to the '-n' option\n"
2326 static gmx_bool do_numbering(int natoms, gmx_groups_t *groups, int ng, char *ptrs[],
2327 t_blocka *block, char *gnames[],
2328 int gtype, int restnm,
2329 int grptp, gmx_bool bVerbose,
2332 unsigned short *cbuf;
2333 t_grps *grps = &(groups->grps[gtype]);
2334 int i, j, gid, aj, ognr, ntot = 0;
2337 char warn_buf[STRLEN];
2341 fprintf(debug, "Starting numbering %d groups of type %d\n", ng, gtype);
2344 title = gtypes[gtype];
2347 /* Mark all id's as not set */
2348 for (i = 0; (i < natoms); i++)
2353 snew(grps->nm_ind, ng+1); /* +1 for possible rest group */
2354 for (i = 0; (i < ng); i++)
2356 /* Lookup the group name in the block structure */
2357 gid = search_string(ptrs[i], block->nr, gnames);
2358 if ((grptp != egrptpONE) || (i == 0))
2360 grps->nm_ind[grps->nr++] = gid;
2364 fprintf(debug, "Found gid %d for group %s\n", gid, ptrs[i]);
2367 /* Now go over the atoms in the group */
2368 for (j = block->index[gid]; (j < block->index[gid+1]); j++)
2373 /* Range checking */
2374 if ((aj < 0) || (aj >= natoms))
2376 gmx_fatal(FARGS, "Invalid atom number %d in indexfile", aj);
2378 /* Lookup up the old group number */
2382 gmx_fatal(FARGS, "Atom %d in multiple %s groups (%d and %d)",
2383 aj+1, title, ognr+1, i+1);
2387 /* Store the group number in buffer */
2388 if (grptp == egrptpONE)
2401 /* Now check whether we have done all atoms */
2405 if (grptp == egrptpALL)
2407 gmx_fatal(FARGS, "%d atoms are not part of any of the %s groups",
2408 natoms-ntot, title);
2410 else if (grptp == egrptpPART)
2412 sprintf(warn_buf, "%d atoms are not part of any of the %s groups",
2413 natoms-ntot, title);
2414 warning_note(wi, warn_buf);
2416 /* Assign all atoms currently unassigned to a rest group */
2417 for (j = 0; (j < natoms); j++)
2419 if (cbuf[j] == NOGID)
2425 if (grptp != egrptpPART)
2430 "Making dummy/rest group for %s containing %d elements\n",
2431 title, natoms-ntot);
2433 /* Add group name "rest" */
2434 grps->nm_ind[grps->nr] = restnm;
2436 /* Assign the rest name to all atoms not currently assigned to a group */
2437 for (j = 0; (j < natoms); j++)
2439 if (cbuf[j] == NOGID)
2448 if (grps->nr == 1 && (ntot == 0 || ntot == natoms))
2450 /* All atoms are part of one (or no) group, no index required */
2451 groups->ngrpnr[gtype] = 0;
2452 groups->grpnr[gtype] = NULL;
2456 groups->ngrpnr[gtype] = natoms;
2457 snew(groups->grpnr[gtype], natoms);
2458 for (j = 0; (j < natoms); j++)
2460 groups->grpnr[gtype][j] = cbuf[j];
2466 return (bRest && grptp == egrptpPART);
2469 static void calc_nrdf(gmx_mtop_t *mtop, t_inputrec *ir, char **gnames)
2472 gmx_groups_t *groups;
2474 int natoms, ai, aj, i, j, d, g, imin, jmin, nc;
2476 int *nrdf2, *na_vcm, na_tot;
2477 double *nrdf_tc, *nrdf_vcm, nrdf_uc, n_sub = 0;
2478 gmx_mtop_atomloop_all_t aloop;
2480 int mb, mol, ftype, as;
2481 gmx_molblock_t *molb;
2482 gmx_moltype_t *molt;
2485 * First calc 3xnr-atoms for each group
2486 * then subtract half a degree of freedom for each constraint
2488 * Only atoms and nuclei contribute to the degrees of freedom...
2493 groups = &mtop->groups;
2494 natoms = mtop->natoms;
2496 /* Allocate one more for a possible rest group */
2497 /* We need to sum degrees of freedom into doubles,
2498 * since floats give too low nrdf's above 3 million atoms.
2500 snew(nrdf_tc, groups->grps[egcTC].nr+1);
2501 snew(nrdf_vcm, groups->grps[egcVCM].nr+1);
2502 snew(na_vcm, groups->grps[egcVCM].nr+1);
2504 for (i = 0; i < groups->grps[egcTC].nr; i++)
2508 for (i = 0; i < groups->grps[egcVCM].nr+1; i++)
2513 snew(nrdf2, natoms);
2514 aloop = gmx_mtop_atomloop_all_init(mtop);
2515 while (gmx_mtop_atomloop_all_next(aloop, &i, &atom))
2518 if (atom->ptype == eptAtom || atom->ptype == eptNucleus)
2520 g = ggrpnr(groups, egcFREEZE, i);
2521 /* Double count nrdf for particle i */
2522 for (d = 0; d < DIM; d++)
2524 if (opts->nFreeze[g][d] == 0)
2529 nrdf_tc [ggrpnr(groups, egcTC, i)] += 0.5*nrdf2[i];
2530 nrdf_vcm[ggrpnr(groups, egcVCM, i)] += 0.5*nrdf2[i];
2535 for (mb = 0; mb < mtop->nmolblock; mb++)
2537 molb = &mtop->molblock[mb];
2538 molt = &mtop->moltype[molb->type];
2539 atom = molt->atoms.atom;
2540 for (mol = 0; mol < molb->nmol; mol++)
2542 for (ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++)
2544 ia = molt->ilist[ftype].iatoms;
2545 for (i = 0; i < molt->ilist[ftype].nr; )
2547 /* Subtract degrees of freedom for the constraints,
2548 * if the particles still have degrees of freedom left.
2549 * If one of the particles is a vsite or a shell, then all
2550 * constraint motion will go there, but since they do not
2551 * contribute to the constraints the degrees of freedom do not
2556 if (((atom[ia[1]].ptype == eptNucleus) ||
2557 (atom[ia[1]].ptype == eptAtom)) &&
2558 ((atom[ia[2]].ptype == eptNucleus) ||
2559 (atom[ia[2]].ptype == eptAtom)))
2577 imin = min(imin, nrdf2[ai]);
2578 jmin = min(jmin, nrdf2[aj]);
2581 nrdf_tc [ggrpnr(groups, egcTC, ai)] -= 0.5*imin;
2582 nrdf_tc [ggrpnr(groups, egcTC, aj)] -= 0.5*jmin;
2583 nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5*imin;
2584 nrdf_vcm[ggrpnr(groups, egcVCM, aj)] -= 0.5*jmin;
2586 ia += interaction_function[ftype].nratoms+1;
2587 i += interaction_function[ftype].nratoms+1;
2590 ia = molt->ilist[F_SETTLE].iatoms;
2591 for (i = 0; i < molt->ilist[F_SETTLE].nr; )
2593 /* Subtract 1 dof from every atom in the SETTLE */
2594 for (j = 0; j < 3; j++)
2597 imin = min(2, nrdf2[ai]);
2599 nrdf_tc [ggrpnr(groups, egcTC, ai)] -= 0.5*imin;
2600 nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5*imin;
2605 as += molt->atoms.nr;
2609 if (ir->ePull == epullCONSTRAINT)
2611 /* Correct nrdf for the COM constraints.
2612 * We correct using the TC and VCM group of the first atom
2613 * in the reference and pull group. If atoms in one pull group
2614 * belong to different TC or VCM groups it is anyhow difficult
2615 * to determine the optimal nrdf assignment.
2618 if (pull->eGeom == epullgPOS)
2621 for (i = 0; i < DIM; i++)
2633 for (i = 0; i < pull->ngrp; i++)
2636 if (pull->grp[0].nat > 0)
2638 /* Subtract 1/2 dof from the reference group */
2639 ai = pull->grp[0].ind[0];
2640 if (nrdf_tc[ggrpnr(groups, egcTC, ai)] > 1)
2642 nrdf_tc [ggrpnr(groups, egcTC, ai)] -= 0.5;
2643 nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5;
2647 /* Subtract 1/2 dof from the pulled group */
2648 ai = pull->grp[1+i].ind[0];
2649 nrdf_tc [ggrpnr(groups, egcTC, ai)] -= 0.5*imin;
2650 nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5*imin;
2651 if (nrdf_tc[ggrpnr(groups, egcTC, ai)] < 0)
2653 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)]]);
2658 if (ir->nstcomm != 0)
2660 /* Subtract 3 from the number of degrees of freedom in each vcm group
2661 * when com translation is removed and 6 when rotation is removed
2664 switch (ir->comm_mode)
2667 n_sub = ndof_com(ir);
2674 gmx_incons("Checking comm_mode");
2677 for (i = 0; i < groups->grps[egcTC].nr; i++)
2679 /* Count the number of atoms of TC group i for every VCM group */
2680 for (j = 0; j < groups->grps[egcVCM].nr+1; j++)
2685 for (ai = 0; ai < natoms; ai++)
2687 if (ggrpnr(groups, egcTC, ai) == i)
2689 na_vcm[ggrpnr(groups, egcVCM, ai)]++;
2693 /* Correct for VCM removal according to the fraction of each VCM
2694 * group present in this TC group.
2696 nrdf_uc = nrdf_tc[i];
2699 fprintf(debug, "T-group[%d] nrdf_uc = %g, n_sub = %g\n",
2703 for (j = 0; j < groups->grps[egcVCM].nr+1; j++)
2705 if (nrdf_vcm[j] > n_sub)
2707 nrdf_tc[i] += nrdf_uc*((double)na_vcm[j]/(double)na_tot)*
2708 (nrdf_vcm[j] - n_sub)/nrdf_vcm[j];
2712 fprintf(debug, " nrdf_vcm[%d] = %g, nrdf = %g\n",
2713 j, nrdf_vcm[j], nrdf_tc[i]);
2718 for (i = 0; (i < groups->grps[egcTC].nr); i++)
2720 opts->nrdf[i] = nrdf_tc[i];
2721 if (opts->nrdf[i] < 0)
2726 "Number of degrees of freedom in T-Coupling group %s is %.2f\n",
2727 gnames[groups->grps[egcTC].nm_ind[i]], opts->nrdf[i]);
2736 static void decode_cos(char *s, t_cosines *cosine, gmx_bool bTime)
2739 char format[STRLEN], f1[STRLEN];
2751 sscanf(t, "%d", &(cosine->n));
2758 snew(cosine->a, cosine->n);
2759 snew(cosine->phi, cosine->n);
2761 sprintf(format, "%%*d");
2762 for (i = 0; (i < cosine->n); i++)
2765 strcat(f1, "%lf%lf");
2766 if (sscanf(t, f1, &a, &phi) < 2)
2768 gmx_fatal(FARGS, "Invalid input for electric field shift: '%s'", t);
2771 cosine->phi[i] = phi;
2772 strcat(format, "%*lf%*lf");
2779 static gmx_bool do_egp_flag(t_inputrec *ir, gmx_groups_t *groups,
2780 const char *option, const char *val, int flag)
2782 /* The maximum number of energy group pairs would be MAXPTR*(MAXPTR+1)/2.
2783 * But since this is much larger than STRLEN, such a line can not be parsed.
2784 * The real maximum is the number of names that fit in a string: STRLEN/2.
2786 #define EGP_MAX (STRLEN/2)
2787 int nelem, i, j, k, nr;
2788 char *names[EGP_MAX];
2792 gnames = groups->grpname;
2794 nelem = str_nelem(val, EGP_MAX, names);
2797 gmx_fatal(FARGS, "The number of groups for %s is odd", option);
2799 nr = groups->grps[egcENER].nr;
2801 for (i = 0; i < nelem/2; i++)
2805 gmx_strcasecmp(names[2*i], *(gnames[groups->grps[egcENER].nm_ind[j]])))
2811 gmx_fatal(FARGS, "%s in %s is not an energy group\n",
2812 names[2*i], option);
2816 gmx_strcasecmp(names[2*i+1], *(gnames[groups->grps[egcENER].nm_ind[k]])))
2822 gmx_fatal(FARGS, "%s in %s is not an energy group\n",
2823 names[2*i+1], option);
2825 if ((j < nr) && (k < nr))
2827 ir->opts.egp_flags[nr*j+k] |= flag;
2828 ir->opts.egp_flags[nr*k+j] |= flag;
2836 void do_index(const char* mdparin, const char *ndx,
2839 t_inputrec *ir, rvec *v,
2843 gmx_groups_t *groups;
2847 char warnbuf[STRLEN], **gnames;
2848 int nr, ntcg, ntau_t, nref_t, nacc, nofg, nSA, nSA_points, nSA_time, nSA_temp;
2851 int nacg, nfreeze, nfrdim, nenergy, nvcm, nuser;
2852 char *ptr1[MAXPTR], *ptr2[MAXPTR], *ptr3[MAXPTR];
2853 int i, j, k, restnm;
2855 gmx_bool bExcl, bTable, bSetTCpar, bAnneal, bRest;
2856 int nQMmethod, nQMbasis, nQMcharge, nQMmult, nbSH, nCASorb, nCASelec,
2857 nSAon, nSAoff, nSAsteps, nQMg, nbOPT, nbTS;
2858 char warn_buf[STRLEN];
2862 fprintf(stderr, "processing index file...\n");
2868 snew(grps->index, 1);
2870 atoms_all = gmx_mtop_global_atoms(mtop);
2871 analyse(&atoms_all, grps, &gnames, FALSE, TRUE);
2872 free_t_atoms(&atoms_all, FALSE);
2876 grps = init_index(ndx, &gnames);
2879 groups = &mtop->groups;
2880 natoms = mtop->natoms;
2881 symtab = &mtop->symtab;
2883 snew(groups->grpname, grps->nr+1);
2885 for (i = 0; (i < grps->nr); i++)
2887 groups->grpname[i] = put_symtab(symtab, gnames[i]);
2889 groups->grpname[i] = put_symtab(symtab, "rest");
2891 srenew(gnames, grps->nr+1);
2892 gnames[restnm] = *(groups->grpname[i]);
2893 groups->ngrpname = grps->nr+1;
2895 set_warning_line(wi, mdparin, -1);
2897 ntau_t = str_nelem(tau_t, MAXPTR, ptr1);
2898 nref_t = str_nelem(ref_t, MAXPTR, ptr2);
2899 ntcg = str_nelem(tcgrps, MAXPTR, ptr3);
2900 if ((ntau_t != ntcg) || (nref_t != ntcg))
2902 gmx_fatal(FARGS, "Invalid T coupling input: %d groups, %d ref-t values and "
2903 "%d tau-t values", ntcg, nref_t, ntau_t);
2906 bSetTCpar = (ir->etc || EI_SD(ir->eI) || ir->eI == eiBD || EI_TPI(ir->eI));
2907 do_numbering(natoms, groups, ntcg, ptr3, grps, gnames, egcTC,
2908 restnm, bSetTCpar ? egrptpALL : egrptpALL_GENREST, bVerbose, wi);
2909 nr = groups->grps[egcTC].nr;
2911 snew(ir->opts.nrdf, nr);
2912 snew(ir->opts.tau_t, nr);
2913 snew(ir->opts.ref_t, nr);
2914 if (ir->eI == eiBD && ir->bd_fric == 0)
2916 fprintf(stderr, "bd-fric=0, so tau-t will be used as the inverse friction constant(s)\n");
2923 gmx_fatal(FARGS, "Not enough ref-t and tau-t values!");
2927 for (i = 0; (i < nr); i++)
2929 ir->opts.tau_t[i] = strtod(ptr1[i], NULL);
2930 if ((ir->eI == eiBD || ir->eI == eiSD2) && ir->opts.tau_t[i] <= 0)
2932 sprintf(warn_buf, "With integrator %s tau-t should be larger than 0", ei_names[ir->eI]);
2933 warning_error(wi, warn_buf);
2936 if (ir->etc != etcVRESCALE && ir->opts.tau_t[i] == 0)
2938 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.");
2941 if (ir->opts.tau_t[i] >= 0)
2943 tau_min = min(tau_min, ir->opts.tau_t[i]);
2946 if (ir->etc != etcNO && ir->nsttcouple == -1)
2948 ir->nsttcouple = ir_optimal_nsttcouple(ir);
2953 if ((ir->etc == etcNOSEHOOVER) && (ir->epc == epcBERENDSEN))
2955 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");
2957 if ((ir->epc == epcMTTK) && (ir->etc > etcNO))
2959 if (ir->nstpcouple != ir->nsttcouple)
2961 int mincouple = min(ir->nstpcouple, ir->nsttcouple);
2962 ir->nstpcouple = ir->nsttcouple = mincouple;
2963 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);
2964 warning_note(wi, warn_buf);
2968 /* velocity verlet with averaged kinetic energy KE = 0.5*(v(t+1/2) - v(t-1/2)) is implemented
2969 primarily for testing purposes, and does not work with temperature coupling other than 1 */
2971 if (ETC_ANDERSEN(ir->etc))
2973 if (ir->nsttcouple != 1)
2976 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");
2977 warning_note(wi, warn_buf);
2980 nstcmin = tcouple_min_integration_steps(ir->etc);
2983 if (tau_min/(ir->delta_t*ir->nsttcouple) < nstcmin)
2985 sprintf(warn_buf, "For proper integration of the %s thermostat, tau-t (%g) should be at least %d times larger than nsttcouple*dt (%g)",
2986 ETCOUPLTYPE(ir->etc),
2988 ir->nsttcouple*ir->delta_t);
2989 warning(wi, warn_buf);
2992 for (i = 0; (i < nr); i++)
2994 ir->opts.ref_t[i] = strtod(ptr2[i], NULL);
2995 if (ir->opts.ref_t[i] < 0)
2997 gmx_fatal(FARGS, "ref-t for group %d negative", i);
3000 /* set the lambda mc temperature to the md integrator temperature (which should be defined
3001 if we are in this conditional) if mc_temp is negative */
3002 if (ir->expandedvals->mc_temp < 0)
3004 ir->expandedvals->mc_temp = ir->opts.ref_t[0]; /*for now, set to the first reft */
3008 /* Simulated annealing for each group. There are nr groups */
3009 nSA = str_nelem(anneal, MAXPTR, ptr1);
3010 if (nSA == 1 && (ptr1[0][0] == 'n' || ptr1[0][0] == 'N'))
3014 if (nSA > 0 && nSA != nr)
3016 gmx_fatal(FARGS, "Not enough annealing values: %d (for %d groups)\n", nSA, nr);
3020 snew(ir->opts.annealing, nr);
3021 snew(ir->opts.anneal_npoints, nr);
3022 snew(ir->opts.anneal_time, nr);
3023 snew(ir->opts.anneal_temp, nr);
3024 for (i = 0; i < nr; i++)
3026 ir->opts.annealing[i] = eannNO;
3027 ir->opts.anneal_npoints[i] = 0;
3028 ir->opts.anneal_time[i] = NULL;
3029 ir->opts.anneal_temp[i] = NULL;
3034 for (i = 0; i < nr; i++)
3036 if (ptr1[i][0] == 'n' || ptr1[i][0] == 'N')
3038 ir->opts.annealing[i] = eannNO;
3040 else if (ptr1[i][0] == 's' || ptr1[i][0] == 'S')
3042 ir->opts.annealing[i] = eannSINGLE;
3045 else if (ptr1[i][0] == 'p' || ptr1[i][0] == 'P')
3047 ir->opts.annealing[i] = eannPERIODIC;
3053 /* Read the other fields too */
3054 nSA_points = str_nelem(anneal_npoints, MAXPTR, ptr1);
3055 if (nSA_points != nSA)
3057 gmx_fatal(FARGS, "Found %d annealing-npoints values for %d groups\n", nSA_points, nSA);
3059 for (k = 0, i = 0; i < nr; i++)
3061 ir->opts.anneal_npoints[i] = strtol(ptr1[i], NULL, 10);
3062 if (ir->opts.anneal_npoints[i] == 1)
3064 gmx_fatal(FARGS, "Please specify at least a start and an end point for annealing\n");
3066 snew(ir->opts.anneal_time[i], ir->opts.anneal_npoints[i]);
3067 snew(ir->opts.anneal_temp[i], ir->opts.anneal_npoints[i]);
3068 k += ir->opts.anneal_npoints[i];
3071 nSA_time = str_nelem(anneal_time, MAXPTR, ptr1);
3074 gmx_fatal(FARGS, "Found %d annealing-time values, wanter %d\n", nSA_time, k);
3076 nSA_temp = str_nelem(anneal_temp, MAXPTR, ptr2);
3079 gmx_fatal(FARGS, "Found %d annealing-temp values, wanted %d\n", nSA_temp, k);
3082 for (i = 0, k = 0; i < nr; i++)
3085 for (j = 0; j < ir->opts.anneal_npoints[i]; j++)
3087 ir->opts.anneal_time[i][j] = strtod(ptr1[k], NULL);
3088 ir->opts.anneal_temp[i][j] = strtod(ptr2[k], NULL);
3091 if (ir->opts.anneal_time[i][0] > (ir->init_t+GMX_REAL_EPS))
3093 gmx_fatal(FARGS, "First time point for annealing > init_t.\n");
3099 if (ir->opts.anneal_time[i][j] < ir->opts.anneal_time[i][j-1])
3101 gmx_fatal(FARGS, "Annealing timepoints out of order: t=%f comes after t=%f\n",
3102 ir->opts.anneal_time[i][j], ir->opts.anneal_time[i][j-1]);
3105 if (ir->opts.anneal_temp[i][j] < 0)
3107 gmx_fatal(FARGS, "Found negative temperature in annealing: %f\n", ir->opts.anneal_temp[i][j]);
3112 /* Print out some summary information, to make sure we got it right */
3113 for (i = 0, k = 0; i < nr; i++)
3115 if (ir->opts.annealing[i] != eannNO)
3117 j = groups->grps[egcTC].nm_ind[i];
3118 fprintf(stderr, "Simulated annealing for group %s: %s, %d timepoints\n",
3119 *(groups->grpname[j]), eann_names[ir->opts.annealing[i]],
3120 ir->opts.anneal_npoints[i]);
3121 fprintf(stderr, "Time (ps) Temperature (K)\n");
3122 /* All terms except the last one */
3123 for (j = 0; j < (ir->opts.anneal_npoints[i]-1); j++)
3125 fprintf(stderr, "%9.1f %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3128 /* Finally the last one */
3129 j = ir->opts.anneal_npoints[i]-1;
3130 if (ir->opts.annealing[i] == eannSINGLE)
3132 fprintf(stderr, "%9.1f- %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3136 fprintf(stderr, "%9.1f %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3137 if (fabs(ir->opts.anneal_temp[i][j]-ir->opts.anneal_temp[i][0]) > GMX_REAL_EPS)
3139 warning_note(wi, "There is a temperature jump when your annealing loops back.\n");
3148 if (ir->ePull != epullNO)
3150 make_pull_groups(ir->pull, pull_grp, grps, gnames);
3155 make_rotation_groups(ir->rot, rot_grp, grps, gnames);
3158 nacc = str_nelem(acc, MAXPTR, ptr1);
3159 nacg = str_nelem(accgrps, MAXPTR, ptr2);
3160 if (nacg*DIM != nacc)
3162 gmx_fatal(FARGS, "Invalid Acceleration input: %d groups and %d acc. values",
3165 do_numbering(natoms, groups, nacg, ptr2, grps, gnames, egcACC,
3166 restnm, egrptpALL_GENREST, bVerbose, wi);
3167 nr = groups->grps[egcACC].nr;
3168 snew(ir->opts.acc, nr);
3169 ir->opts.ngacc = nr;
3171 for (i = k = 0; (i < nacg); i++)
3173 for (j = 0; (j < DIM); j++, k++)
3175 ir->opts.acc[i][j] = strtod(ptr1[k], NULL);
3178 for (; (i < nr); i++)
3180 for (j = 0; (j < DIM); j++)
3182 ir->opts.acc[i][j] = 0;
3186 nfrdim = str_nelem(frdim, MAXPTR, ptr1);
3187 nfreeze = str_nelem(freeze, MAXPTR, ptr2);
3188 if (nfrdim != DIM*nfreeze)
3190 gmx_fatal(FARGS, "Invalid Freezing input: %d groups and %d freeze values",
3193 do_numbering(natoms, groups, nfreeze, ptr2, grps, gnames, egcFREEZE,
3194 restnm, egrptpALL_GENREST, bVerbose, wi);
3195 nr = groups->grps[egcFREEZE].nr;
3196 ir->opts.ngfrz = nr;
3197 snew(ir->opts.nFreeze, nr);
3198 for (i = k = 0; (i < nfreeze); i++)
3200 for (j = 0; (j < DIM); j++, k++)
3202 ir->opts.nFreeze[i][j] = (gmx_strncasecmp(ptr1[k], "Y", 1) == 0);
3203 if (!ir->opts.nFreeze[i][j])
3205 if (gmx_strncasecmp(ptr1[k], "N", 1) != 0)
3207 sprintf(warnbuf, "Please use Y(ES) or N(O) for freezedim only "
3208 "(not %s)", ptr1[k]);
3209 warning(wi, warn_buf);
3214 for (; (i < nr); i++)
3216 for (j = 0; (j < DIM); j++)
3218 ir->opts.nFreeze[i][j] = 0;
3222 nenergy = str_nelem(energy, MAXPTR, ptr1);
3223 do_numbering(natoms, groups, nenergy, ptr1, grps, gnames, egcENER,
3224 restnm, egrptpALL_GENREST, bVerbose, wi);
3225 add_wall_energrps(groups, ir->nwall, symtab);
3226 ir->opts.ngener = groups->grps[egcENER].nr;
3227 nvcm = str_nelem(vcm, MAXPTR, ptr1);
3229 do_numbering(natoms, groups, nvcm, ptr1, grps, gnames, egcVCM,
3230 restnm, nvcm == 0 ? egrptpALL_GENREST : egrptpPART, bVerbose, wi);
3233 warning(wi, "Some atoms are not part of any center of mass motion removal group.\n"
3234 "This may lead to artifacts.\n"
3235 "In most cases one should use one group for the whole system.");
3238 /* Now we have filled the freeze struct, so we can calculate NRDF */
3239 calc_nrdf(mtop, ir, gnames);
3245 /* Must check per group! */
3246 for (i = 0; (i < ir->opts.ngtc); i++)
3248 ntot += ir->opts.nrdf[i];
3250 if (ntot != (DIM*natoms))
3252 fac = sqrt(ntot/(DIM*natoms));
3255 fprintf(stderr, "Scaling velocities by a factor of %.3f to account for constraints\n"
3256 "and removal of center of mass motion\n", fac);
3258 for (i = 0; (i < natoms); i++)
3260 svmul(fac, v[i], v[i]);
3265 nuser = str_nelem(user1, MAXPTR, ptr1);
3266 do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcUser1,
3267 restnm, egrptpALL_GENREST, bVerbose, wi);
3268 nuser = str_nelem(user2, MAXPTR, ptr1);
3269 do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcUser2,
3270 restnm, egrptpALL_GENREST, bVerbose, wi);
3271 nuser = str_nelem(xtc_grps, MAXPTR, ptr1);
3272 do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcXTC,
3273 restnm, egrptpONE, bVerbose, wi);
3274 nofg = str_nelem(orirefitgrp, MAXPTR, ptr1);
3275 do_numbering(natoms, groups, nofg, ptr1, grps, gnames, egcORFIT,
3276 restnm, egrptpALL_GENREST, bVerbose, wi);
3278 /* QMMM input processing */
3279 nQMg = str_nelem(QMMM, MAXPTR, ptr1);
3280 nQMmethod = str_nelem(QMmethod, MAXPTR, ptr2);
3281 nQMbasis = str_nelem(QMbasis, MAXPTR, ptr3);
3282 if ((nQMmethod != nQMg) || (nQMbasis != nQMg))
3284 gmx_fatal(FARGS, "Invalid QMMM input: %d groups %d basissets"
3285 " and %d methods\n", nQMg, nQMbasis, nQMmethod);
3287 /* group rest, if any, is always MM! */
3288 do_numbering(natoms, groups, nQMg, ptr1, grps, gnames, egcQMMM,
3289 restnm, egrptpALL_GENREST, bVerbose, wi);
3290 nr = nQMg; /*atoms->grps[egcQMMM].nr;*/
3291 ir->opts.ngQM = nQMg;
3292 snew(ir->opts.QMmethod, nr);
3293 snew(ir->opts.QMbasis, nr);
3294 for (i = 0; i < nr; i++)
3296 /* input consists of strings: RHF CASSCF PM3 .. These need to be
3297 * converted to the corresponding enum in names.c
3299 ir->opts.QMmethod[i] = search_QMstring(ptr2[i], eQMmethodNR,
3301 ir->opts.QMbasis[i] = search_QMstring(ptr3[i], eQMbasisNR,
3305 nQMmult = str_nelem(QMmult, MAXPTR, ptr1);
3306 nQMcharge = str_nelem(QMcharge, MAXPTR, ptr2);
3307 nbSH = str_nelem(bSH, MAXPTR, ptr3);
3308 snew(ir->opts.QMmult, nr);
3309 snew(ir->opts.QMcharge, nr);
3310 snew(ir->opts.bSH, nr);
3312 for (i = 0; i < nr; i++)
3314 ir->opts.QMmult[i] = strtol(ptr1[i], NULL, 10);
3315 ir->opts.QMcharge[i] = strtol(ptr2[i], NULL, 10);
3316 ir->opts.bSH[i] = (gmx_strncasecmp(ptr3[i], "Y", 1) == 0);
3319 nCASelec = str_nelem(CASelectrons, MAXPTR, ptr1);
3320 nCASorb = str_nelem(CASorbitals, MAXPTR, ptr2);
3321 snew(ir->opts.CASelectrons, nr);
3322 snew(ir->opts.CASorbitals, nr);
3323 for (i = 0; i < nr; i++)
3325 ir->opts.CASelectrons[i] = strtol(ptr1[i], NULL, 10);
3326 ir->opts.CASorbitals[i] = strtol(ptr2[i], NULL, 10);
3328 /* special optimization options */
3330 nbOPT = str_nelem(bOPT, MAXPTR, ptr1);
3331 nbTS = str_nelem(bTS, MAXPTR, ptr2);
3332 snew(ir->opts.bOPT, nr);
3333 snew(ir->opts.bTS, nr);
3334 for (i = 0; i < nr; i++)
3336 ir->opts.bOPT[i] = (gmx_strncasecmp(ptr1[i], "Y", 1) == 0);
3337 ir->opts.bTS[i] = (gmx_strncasecmp(ptr2[i], "Y", 1) == 0);
3339 nSAon = str_nelem(SAon, MAXPTR, ptr1);
3340 nSAoff = str_nelem(SAoff, MAXPTR, ptr2);
3341 nSAsteps = str_nelem(SAsteps, MAXPTR, ptr3);
3342 snew(ir->opts.SAon, nr);
3343 snew(ir->opts.SAoff, nr);
3344 snew(ir->opts.SAsteps, nr);
3346 for (i = 0; i < nr; i++)
3348 ir->opts.SAon[i] = strtod(ptr1[i], NULL);
3349 ir->opts.SAoff[i] = strtod(ptr2[i], NULL);
3350 ir->opts.SAsteps[i] = strtol(ptr3[i], NULL, 10);
3352 /* end of QMMM input */
3356 for (i = 0; (i < egcNR); i++)
3358 fprintf(stderr, "%-16s has %d element(s):", gtypes[i], groups->grps[i].nr);
3359 for (j = 0; (j < groups->grps[i].nr); j++)
3361 fprintf(stderr, " %s", *(groups->grpname[groups->grps[i].nm_ind[j]]));
3363 fprintf(stderr, "\n");
3367 nr = groups->grps[egcENER].nr;
3368 snew(ir->opts.egp_flags, nr*nr);
3370 bExcl = do_egp_flag(ir, groups, "energygrp-excl", egpexcl, EGP_EXCL);
3371 if (bExcl && ir->cutoff_scheme == ecutsVERLET)
3373 warning_error(wi, "Energy group exclusions are not (yet) implemented for the Verlet scheme");
3375 if (bExcl && EEL_FULL(ir->coulombtype))
3377 warning(wi, "Can not exclude the lattice Coulomb energy between energy groups");
3380 bTable = do_egp_flag(ir, groups, "energygrp-table", egptable, EGP_TABLE);
3381 if (bTable && !(ir->vdwtype == evdwUSER) &&
3382 !(ir->coulombtype == eelUSER) && !(ir->coulombtype == eelPMEUSER) &&
3383 !(ir->coulombtype == eelPMEUSERSWITCH))
3385 gmx_fatal(FARGS, "Can only have energy group pair tables in combination with user tables for VdW and/or Coulomb");
3388 decode_cos(efield_x, &(ir->ex[XX]), FALSE);
3389 decode_cos(efield_xt, &(ir->et[XX]), TRUE);
3390 decode_cos(efield_y, &(ir->ex[YY]), FALSE);
3391 decode_cos(efield_yt, &(ir->et[YY]), TRUE);
3392 decode_cos(efield_z, &(ir->ex[ZZ]), FALSE);
3393 decode_cos(efield_zt, &(ir->et[ZZ]), TRUE);
3397 do_adress_index(ir->adress, groups, gnames, &(ir->opts), wi);
3400 for (i = 0; (i < grps->nr); i++)
3412 static void check_disre(gmx_mtop_t *mtop)
3414 gmx_ffparams_t *ffparams;
3415 t_functype *functype;
3417 int i, ndouble, ftype;
3418 int label, old_label;
3420 if (gmx_mtop_ftype_count(mtop, F_DISRES) > 0)
3422 ffparams = &mtop->ffparams;
3423 functype = ffparams->functype;
3424 ip = ffparams->iparams;
3427 for (i = 0; i < ffparams->ntypes; i++)
3429 ftype = functype[i];
3430 if (ftype == F_DISRES)
3432 label = ip[i].disres.label;
3433 if (label == old_label)
3435 fprintf(stderr, "Distance restraint index %d occurs twice\n", label);
3443 gmx_fatal(FARGS, "Found %d double distance restraint indices,\n"
3444 "probably the parameters for multiple pairs in one restraint "
3445 "are not identical\n", ndouble);
3450 static gmx_bool absolute_reference(t_inputrec *ir, gmx_mtop_t *sys,
3451 gmx_bool posres_only,
3455 gmx_mtop_ilistloop_t iloop;
3465 for (d = 0; d < DIM; d++)
3467 AbsRef[d] = (d < ndof_com(ir) ? 0 : 1);
3469 /* Check for freeze groups */
3470 for (g = 0; g < ir->opts.ngfrz; g++)
3472 for (d = 0; d < DIM; d++)
3474 if (ir->opts.nFreeze[g][d] != 0)
3482 /* Check for position restraints */
3483 iloop = gmx_mtop_ilistloop_init(sys);
3484 while (gmx_mtop_ilistloop_next(iloop, &ilist, &nmol))
3487 (AbsRef[XX] == 0 || AbsRef[YY] == 0 || AbsRef[ZZ] == 0))
3489 for (i = 0; i < ilist[F_POSRES].nr; i += 2)
3491 pr = &sys->ffparams.iparams[ilist[F_POSRES].iatoms[i]];
3492 for (d = 0; d < DIM; d++)
3494 if (pr->posres.fcA[d] != 0)
3500 for (i = 0; i < ilist[F_FBPOSRES].nr; i += 2)
3502 /* Check for flat-bottom posres */
3503 pr = &sys->ffparams.iparams[ilist[F_FBPOSRES].iatoms[i]];
3504 if (pr->fbposres.k != 0)
3506 switch (pr->fbposres.geom)
3508 case efbposresSPHERE:
3509 AbsRef[XX] = AbsRef[YY] = AbsRef[ZZ] = 1;
3511 case efbposresCYLINDER:
3512 AbsRef[XX] = AbsRef[YY] = 1;
3514 case efbposresX: /* d=XX */
3515 case efbposresY: /* d=YY */
3516 case efbposresZ: /* d=ZZ */
3517 d = pr->fbposres.geom - efbposresX;
3521 gmx_fatal(FARGS, " Invalid geometry for flat-bottom position restraint.\n"
3522 "Expected nr between 1 and %d. Found %d\n", efbposresNR-1,
3530 return (AbsRef[XX] != 0 && AbsRef[YY] != 0 && AbsRef[ZZ] != 0);
3533 void triple_check(const char *mdparin, t_inputrec *ir, gmx_mtop_t *sys,
3537 int i, m, g, nmol, npct;
3538 gmx_bool bCharge, bAcc;
3539 real gdt_max, *mgrp, mt;
3541 gmx_mtop_atomloop_block_t aloopb;
3542 gmx_mtop_atomloop_all_t aloop;
3545 char warn_buf[STRLEN];
3547 set_warning_line(wi, mdparin, -1);
3549 if (EI_DYNAMICS(ir->eI) && !EI_SD(ir->eI) && ir->eI != eiBD &&
3550 ir->comm_mode == ecmNO &&
3551 !(absolute_reference(ir, sys, FALSE, AbsRef) || ir->nsteps <= 10))
3553 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");
3556 /* Check for pressure coupling with absolute position restraints */
3557 if (ir->epc != epcNO && ir->refcoord_scaling == erscNO)
3559 absolute_reference(ir, sys, TRUE, AbsRef);
3561 for (m = 0; m < DIM; m++)
3563 if (AbsRef[m] && norm2(ir->compress[m]) > 0)
3565 warning(wi, "You are using pressure coupling with absolute position restraints, this will give artifacts. Use the refcoord_scaling option.");
3573 aloopb = gmx_mtop_atomloop_block_init(sys);
3574 while (gmx_mtop_atomloop_block_next(aloopb, &atom, &nmol))
3576 if (atom->q != 0 || atom->qB != 0)
3584 if (EEL_FULL(ir->coulombtype))
3587 "You are using full electrostatics treatment %s for a system without charges.\n"
3588 "This costs a lot of performance for just processing zeros, consider using %s instead.\n",
3589 EELTYPE(ir->coulombtype), EELTYPE(eelCUT));
3590 warning(wi, err_buf);
3595 if (ir->coulombtype == eelCUT && ir->rcoulomb > 0 && !ir->implicit_solvent)
3598 "You are using a plain Coulomb cut-off, which might produce artifacts.\n"
3599 "You might want to consider using %s electrostatics.\n",
3601 warning_note(wi, err_buf);
3605 /* Generalized reaction field */
3606 if (ir->opts.ngtc == 0)
3608 sprintf(err_buf, "No temperature coupling while using coulombtype %s",
3610 CHECK(ir->coulombtype == eelGRF);
3614 sprintf(err_buf, "When using coulombtype = %s"
3615 " ref-t for temperature coupling should be > 0",
3617 CHECK((ir->coulombtype == eelGRF) && (ir->opts.ref_t[0] <= 0));
3620 if (ir->eI == eiSD1 &&
3621 (gmx_mtop_ftype_count(sys, F_CONSTR) > 0 ||
3622 gmx_mtop_ftype_count(sys, F_SETTLE) > 0))
3624 sprintf(warn_buf, "With constraints integrator %s is less accurate, consider using %s instead", ei_names[ir->eI], ei_names[eiSD2]);
3625 warning_note(wi, warn_buf);
3629 for (i = 0; (i < sys->groups.grps[egcACC].nr); i++)
3631 for (m = 0; (m < DIM); m++)
3633 if (fabs(ir->opts.acc[i][m]) > 1e-6)
3642 snew(mgrp, sys->groups.grps[egcACC].nr);
3643 aloop = gmx_mtop_atomloop_all_init(sys);
3644 while (gmx_mtop_atomloop_all_next(aloop, &i, &atom))
3646 mgrp[ggrpnr(&sys->groups, egcACC, i)] += atom->m;
3649 for (i = 0; (i < sys->groups.grps[egcACC].nr); i++)
3651 for (m = 0; (m < DIM); m++)
3653 acc[m] += ir->opts.acc[i][m]*mgrp[i];
3657 for (m = 0; (m < DIM); m++)
3659 if (fabs(acc[m]) > 1e-6)
3661 const char *dim[DIM] = { "X", "Y", "Z" };
3663 "Net Acceleration in %s direction, will %s be corrected\n",
3664 dim[m], ir->nstcomm != 0 ? "" : "not");
3665 if (ir->nstcomm != 0 && m < ndof_com(ir))
3668 for (i = 0; (i < sys->groups.grps[egcACC].nr); i++)
3670 ir->opts.acc[i][m] -= acc[m];
3678 if (ir->efep != efepNO && ir->fepvals->sc_alpha != 0 &&
3679 !gmx_within_tol(sys->ffparams.reppow, 12.0, 10*GMX_DOUBLE_EPS))
3681 gmx_fatal(FARGS, "Soft-core interactions are only supported with VdW repulsion power 12");
3684 if (ir->ePull != epullNO)
3686 if (ir->pull->grp[0].nat == 0)
3688 absolute_reference(ir, sys, FALSE, AbsRef);
3689 for (m = 0; m < DIM; m++)
3691 if (ir->pull->dim[m] && !AbsRef[m])
3693 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.");
3699 if (ir->pull->eGeom == epullgDIRPBC)
3701 for (i = 0; i < 3; i++)
3703 for (m = 0; m <= i; m++)
3705 if ((ir->epc != epcNO && ir->compress[i][m] != 0) ||
3706 ir->deform[i][m] != 0)
3708 for (g = 1; g < ir->pull->ngrp; g++)
3710 if (ir->pull->grp[g].vec[m] != 0)
3712 gmx_fatal(FARGS, "Can not have dynamic box while using pull geometry '%s' (dim %c)", EPULLGEOM(ir->pull->eGeom), 'x'+m);
3724 void double_check(t_inputrec *ir, matrix box, gmx_bool bConstr, warninp_t wi)
3728 char warn_buf[STRLEN];
3731 ptr = check_box(ir->ePBC, box);
3734 warning_error(wi, ptr);
3737 if (bConstr && ir->eConstrAlg == econtSHAKE)
3739 if (ir->shake_tol <= 0.0)
3741 sprintf(warn_buf, "ERROR: shake-tol must be > 0 instead of %g\n",
3743 warning_error(wi, warn_buf);
3746 if (IR_TWINRANGE(*ir) && ir->nstlist > 1)
3748 sprintf(warn_buf, "With twin-range cut-off's and SHAKE the virial and the pressure are incorrect.");
3749 if (ir->epc == epcNO)
3751 warning(wi, warn_buf);
3755 warning_error(wi, warn_buf);
3760 if ( (ir->eConstrAlg == econtLINCS) && bConstr)
3762 /* If we have Lincs constraints: */
3763 if (ir->eI == eiMD && ir->etc == etcNO &&
3764 ir->eConstrAlg == econtLINCS && ir->nLincsIter == 1)
3766 sprintf(warn_buf, "For energy conservation with LINCS, lincs_iter should be 2 or larger.\n");
3767 warning_note(wi, warn_buf);
3770 if ((ir->eI == eiCG || ir->eI == eiLBFGS) && (ir->nProjOrder < 8))
3772 sprintf(warn_buf, "For accurate %s with LINCS constraints, lincs-order should be 8 or more.", ei_names[ir->eI]);
3773 warning_note(wi, warn_buf);
3775 if (ir->epc == epcMTTK)
3777 warning_error(wi, "MTTK not compatible with lincs -- use shake instead.");
3781 if (ir->LincsWarnAngle > 90.0)
3783 sprintf(warn_buf, "lincs-warnangle can not be larger than 90 degrees, setting it to 90.\n");
3784 warning(wi, warn_buf);
3785 ir->LincsWarnAngle = 90.0;
3788 if (ir->ePBC != epbcNONE)
3790 if (ir->nstlist == 0)
3792 warning(wi, "With nstlist=0 atoms are only put into the box at step 0, therefore drifting atoms might cause the simulation to crash.");
3794 bTWIN = (ir->rlistlong > ir->rlist);
3795 if (ir->ns_type == ensGRID)
3797 if (sqr(ir->rlistlong) >= max_cutoff2(ir->ePBC, box))
3799 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",
3800 bTWIN ? (ir->rcoulomb == ir->rlistlong ? "rcoulomb" : "rvdw") : "rlist");
3801 warning_error(wi, warn_buf);
3806 min_size = min(box[XX][XX], min(box[YY][YY], box[ZZ][ZZ]));
3807 if (2*ir->rlistlong >= min_size)
3809 sprintf(warn_buf, "ERROR: One of the box lengths is smaller than twice the cut-off length. Increase the box size or decrease rlist.");
3810 warning_error(wi, warn_buf);
3813 fprintf(stderr, "Grid search might allow larger cut-off's than simple search with triclinic boxes.");
3820 void check_chargegroup_radii(const gmx_mtop_t *mtop, const t_inputrec *ir,
3824 real rvdw1, rvdw2, rcoul1, rcoul2;
3825 char warn_buf[STRLEN];
3827 calc_chargegroup_radii(mtop, x, &rvdw1, &rvdw2, &rcoul1, &rcoul2);
3831 printf("Largest charge group radii for Van der Waals: %5.3f, %5.3f nm\n",
3836 printf("Largest charge group radii for Coulomb: %5.3f, %5.3f nm\n",
3842 if (rvdw1 + rvdw2 > ir->rlist ||
3843 rcoul1 + rcoul2 > ir->rlist)
3846 "The sum of the two largest charge group radii (%f) "
3847 "is larger than rlist (%f)\n",
3848 max(rvdw1+rvdw2, rcoul1+rcoul2), ir->rlist);
3849 warning(wi, warn_buf);
3853 /* Here we do not use the zero at cut-off macro,
3854 * since user defined interactions might purposely
3855 * not be zero at the cut-off.
3857 if ((EVDW_IS_ZERO_AT_CUTOFF(ir->vdwtype) ||
3858 ir->vdw_modifier != eintmodNONE) &&
3859 rvdw1 + rvdw2 > ir->rlistlong - ir->rvdw)
3861 sprintf(warn_buf, "The sum of the two largest charge group "
3862 "radii (%f) is larger than %s (%f) - rvdw (%f).\n"
3863 "With exact cut-offs, better performance can be "
3864 "obtained with cutoff-scheme = %s, because it "
3865 "does not use charge groups at all.",
3867 ir->rlistlong > ir->rlist ? "rlistlong" : "rlist",
3868 ir->rlistlong, ir->rvdw,
3869 ecutscheme_names[ecutsVERLET]);
3872 warning(wi, warn_buf);
3876 warning_note(wi, warn_buf);
3879 if ((EEL_IS_ZERO_AT_CUTOFF(ir->coulombtype) ||
3880 ir->coulomb_modifier != eintmodNONE) &&
3881 rcoul1 + rcoul2 > ir->rlistlong - ir->rcoulomb)
3883 sprintf(warn_buf, "The sum of the two largest charge group radii (%f) is larger than %s (%f) - rcoulomb (%f).\n"
3884 "With exact cut-offs, better performance can be obtained with cutoff-scheme = %s, because it does not use charge groups at all.",
3886 ir->rlistlong > ir->rlist ? "rlistlong" : "rlist",
3887 ir->rlistlong, ir->rcoulomb,
3888 ecutscheme_names[ecutsVERLET]);
3891 warning(wi, warn_buf);
3895 warning_note(wi, warn_buf);