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"
68 /* Resource parameters
69 * Do not change any of these until you read the instruction
70 * in readinp.h. Some cpp's do not take spaces after the backslash
71 * (like the c-shell), which will give you a very weird compiler
75 static char tcgrps[STRLEN],tau_t[STRLEN],ref_t[STRLEN],
76 acc[STRLEN],accgrps[STRLEN],freeze[STRLEN],frdim[STRLEN],
77 energy[STRLEN],user1[STRLEN],user2[STRLEN],vcm[STRLEN],xtc_grps[STRLEN],
78 couple_moltype[STRLEN],orirefitgrp[STRLEN],egptable[STRLEN],egpexcl[STRLEN],
79 wall_atomtype[STRLEN],wall_density[STRLEN],deform[STRLEN],QMMM[STRLEN];
80 static char foreign_lambda[STRLEN];
81 static char **pull_grp;
82 static char anneal[STRLEN],anneal_npoints[STRLEN],
83 anneal_time[STRLEN],anneal_temp[STRLEN];
84 static char QMmethod[STRLEN],QMbasis[STRLEN],QMcharge[STRLEN],QMmult[STRLEN],
85 bSH[STRLEN],CASorbitals[STRLEN], CASelectrons[STRLEN],SAon[STRLEN],
86 SAoff[STRLEN],SAsteps[STRLEN],bTS[STRLEN],bOPT[STRLEN];
87 static char efield_x[STRLEN],efield_xt[STRLEN],efield_y[STRLEN],
88 efield_yt[STRLEN],efield_z[STRLEN],efield_zt[STRLEN];
91 egrptpALL, /* All particles have to be a member of a group. */
92 egrptpALL_GENREST, /* A rest group with name is generated for particles *
93 * that are not part of any group. */
94 egrptpPART, /* As egrptpALL_GENREST, but no name is generated *
95 * for the rest group. */
96 egrptpONE /* Merge all selected groups into one group, *
97 * make a rest group for the remaining particles. */
101 void init_ir(t_inputrec *ir, t_gromppopts *opts)
103 snew(opts->include,STRLEN);
104 snew(opts->define,STRLEN);
107 static void _low_check(gmx_bool b,char *s,warninp_t wi)
115 static void check_nst(const char *desc_nst,int nst,
116 const char *desc_p,int *p,
121 if (*p > 0 && *p % nst != 0)
123 /* Round up to the next multiple of nst */
124 *p = ((*p)/nst + 1)*nst;
125 sprintf(buf,"%s should be a multiple of %s, changing %s to %d\n",
126 desc_p,desc_nst,desc_p,*p);
131 static gmx_bool ir_NVE(const t_inputrec *ir)
133 return ((ir->eI == eiMD || EI_VV(ir->eI)) && ir->etc == etcNO);
136 static int lcd(int n1,int n2)
141 for(i=2; (i<=n1 && i<=n2); i++)
143 if (n1 % i == 0 && n2 % i == 0)
152 void check_ir(const char *mdparin,t_inputrec *ir, t_gromppopts *opts,
154 /* Check internal consistency */
156 /* Strange macro: first one fills the err_buf, and then one can check
157 * the condition, which will print the message and increase the error
160 #define CHECK(b) _low_check(b,err_buf,wi)
161 char err_buf[256],warn_buf[STRLEN];
165 set_warning_line(wi,mdparin,-1);
167 /* BASIC CUT-OFF STUFF */
168 if (ir->rlist == 0 ||
169 !((EEL_MIGHT_BE_ZERO_AT_CUTOFF(ir->coulombtype) && ir->rcoulomb > ir->rlist) ||
170 (EVDW_MIGHT_BE_ZERO_AT_CUTOFF(ir->vdwtype) && ir->rvdw > ir->rlist))) {
171 /* No switched potential and/or no twin-range:
172 * we can set the long-range cut-off to the maximum of the other cut-offs.
174 ir->rlistlong = max_cutoff(ir->rlist,max_cutoff(ir->rvdw,ir->rcoulomb));
175 } else if (ir->rlistlong < 0) {
176 ir->rlistlong = max_cutoff(ir->rlist,max_cutoff(ir->rvdw,ir->rcoulomb));
177 sprintf(warn_buf,"rlistlong was not set, setting it to %g (no buffer)",
179 warning(wi,warn_buf);
181 if (ir->rlistlong == 0 && ir->ePBC != epbcNONE) {
182 warning_error(wi,"Can not have an infinite cut-off with PBC");
184 if (ir->rlistlong > 0 && (ir->rlist == 0 || ir->rlistlong < ir->rlist)) {
185 warning_error(wi,"rlistlong can not be shorter than rlist");
187 if (IR_TWINRANGE(*ir) && ir->nstlist <= 0) {
188 warning_error(wi,"Can not have nstlist<=0 with twin-range interactions");
191 /* GENERAL INTEGRATOR STUFF */
192 if (!(ir->eI == eiMD || EI_VV(ir->eI)))
196 if (!EI_DYNAMICS(ir->eI))
200 if (EI_DYNAMICS(ir->eI))
202 if (ir->nstcalcenergy < 0)
204 ir->nstcalcenergy = ir_optimal_nstcalcenergy(ir);
205 if (ir->nstenergy != 0 && ir->nstenergy < ir->nstcalcenergy)
207 /* nstcalcenergy larger than nstener does not make sense.
208 * We ideally want nstcalcenergy=nstener.
212 ir->nstcalcenergy = lcd(ir->nstenergy,ir->nstlist);
216 ir->nstcalcenergy = ir->nstenergy;
220 if (ir->epc != epcNO)
222 if (ir->nstpcouple < 0)
224 ir->nstpcouple = ir_optimal_nstpcouple(ir);
227 if (IR_TWINRANGE(*ir))
229 check_nst("nstlist",ir->nstlist,
230 "nstcalcenergy",&ir->nstcalcenergy,wi);
231 if (ir->epc != epcNO)
233 check_nst("nstlist",ir->nstlist,
234 "nstpcouple",&ir->nstpcouple,wi);
238 if (ir->nstcalcenergy > 1)
240 /* for storing exact averages nstenergy should be
241 * a multiple of nstcalcenergy
243 check_nst("nstcalcenergy",ir->nstcalcenergy,
244 "nstenergy",&ir->nstenergy,wi);
245 if (ir->efep != efepNO)
247 /* nstdhdl should be a multiple of nstcalcenergy */
248 check_nst("nstcalcenergy",ir->nstcalcenergy,
249 "nstdhdl",&ir->nstdhdl,wi);
255 if ((EI_SD(ir->eI) || ir->eI == eiBD) &&
256 ir->bContinuation && ir->ld_seed != -1) {
257 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)");
261 if (EI_TPI(ir->eI)) {
262 sprintf(err_buf,"TPI only works with pbc = %s",epbc_names[epbcXYZ]);
263 CHECK(ir->ePBC != epbcXYZ);
264 sprintf(err_buf,"TPI only works with ns = %s",ens_names[ensGRID]);
265 CHECK(ir->ns_type != ensGRID);
266 sprintf(err_buf,"with TPI nstlist should be larger than zero");
267 CHECK(ir->nstlist <= 0);
268 sprintf(err_buf,"TPI does not work with full electrostatics other than PME");
269 CHECK(EEL_FULL(ir->coulombtype) && !EEL_PME(ir->coulombtype));
273 if ( (opts->nshake > 0) && (opts->bMorse) ) {
275 "Using morse bond-potentials while constraining bonds is useless");
276 warning(wi,warn_buf);
279 sprintf(err_buf,"shake_tol must be > 0 instead of %g while using shake",
281 CHECK(((ir->shake_tol <= 0.0) && (opts->nshake>0) &&
282 (ir->eConstrAlg == econtSHAKE)));
285 sprintf(err_buf,"walls only work with pbc=%s",epbc_names[epbcXY]);
286 CHECK(ir->nwall && ir->ePBC!=epbcXY);
289 if (ir->ePBC != epbcXYZ && ir->nwall != 2) {
290 if (ir->ePBC == epbcNONE) {
291 if (ir->epc != epcNO) {
292 warning(wi,"Turning off pressure coupling for vacuum system");
296 sprintf(err_buf,"Can not have pressure coupling with pbc=%s",
297 epbc_names[ir->ePBC]);
298 CHECK(ir->epc != epcNO);
300 sprintf(err_buf,"Can not have Ewald with pbc=%s",epbc_names[ir->ePBC]);
301 CHECK(EEL_FULL(ir->coulombtype));
303 sprintf(err_buf,"Can not have dispersion correction with pbc=%s",
304 epbc_names[ir->ePBC]);
305 CHECK(ir->eDispCorr != edispcNO);
308 if (ir->rlist == 0.0) {
309 sprintf(err_buf,"can only have neighborlist cut-off zero (=infinite)\n"
310 "with coulombtype = %s or coulombtype = %s\n"
311 "without periodic boundary conditions (pbc = %s) and\n"
312 "rcoulomb and rvdw set to zero",
313 eel_names[eelCUT],eel_names[eelUSER],epbc_names[epbcNONE]);
314 CHECK(((ir->coulombtype != eelCUT) && (ir->coulombtype != eelUSER)) ||
315 (ir->ePBC != epbcNONE) ||
316 (ir->rcoulomb != 0.0) || (ir->rvdw != 0.0));
318 if (ir->nstlist < 0) {
319 warning_error(wi,"Can not have heuristic neighborlist updates without cut-off");
321 if (ir->nstlist > 0) {
322 warning_note(wi,"Simulating without cut-offs is usually (slightly) faster with nstlist=0, nstype=simple and particle decomposition");
327 if (ir->nstcomm == 0) {
328 ir->comm_mode = ecmNO;
330 if (ir->comm_mode != ecmNO) {
331 if (ir->nstcomm < 0) {
332 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");
333 ir->nstcomm = abs(ir->nstcomm);
336 if (ir->nstcalcenergy > 0 && ir->nstcomm < ir->nstcalcenergy) {
337 warning_note(wi,"nstcomm < nstcalcenergy defeats the purpose of nstcalcenergy, setting nstcomm to nstcalcenergy");
338 ir->nstcomm = ir->nstcalcenergy;
341 if (ir->comm_mode == ecmANGULAR) {
342 sprintf(err_buf,"Can not remove the rotation around the center of mass with periodic molecules");
343 CHECK(ir->bPeriodicMols);
344 if (ir->ePBC != epbcNONE)
345 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).");
349 if (EI_STATE_VELOCITY(ir->eI) && ir->ePBC == epbcNONE && ir->comm_mode != ecmANGULAR) {
350 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.");
353 sprintf(err_buf,"Free-energy not implemented for Ewald and PPPM");
354 CHECK((ir->coulombtype==eelEWALD || ir->coulombtype==eelPPPM)
355 && (ir->efep!=efepNO));
357 sprintf(err_buf,"Twin-range neighbour searching (NS) with simple NS"
358 " algorithm not implemented");
359 CHECK(((ir->rcoulomb > ir->rlist) || (ir->rvdw > ir->rlist))
360 && (ir->ns_type == ensSIMPLE));
362 /* TEMPERATURE COUPLING */
363 if (ir->etc == etcYES)
365 ir->etc = etcBERENDSEN;
366 warning_note(wi,"Old option for temperature coupling given: "
367 "changing \"yes\" to \"Berendsen\"\n");
370 if (ir->etc == etcNOSEHOOVER)
372 if (ir->opts.nhchainlength < 1)
374 sprintf(warn_buf,"number of Nose-Hoover chains (currently %d) cannot be less than 1,reset to 1\n",ir->opts.nhchainlength);
375 ir->opts.nhchainlength =1;
376 warning(wi,warn_buf);
379 if (ir->etc==etcNOSEHOOVER && !EI_VV(ir->eI) && ir->opts.nhchainlength > 1)
381 warning_note(wi,"leapfrog does not yet support Nose-Hoover chains, nhchainlength reset to 1");
382 ir->opts.nhchainlength = 1;
387 ir->opts.nhchainlength = 0;
390 if (ir->etc == etcBERENDSEN)
392 sprintf(warn_buf,"The %s thermostat does not generate the correct kinetic energy distribution. You might want to consider using the %s thermostat.",
393 ETCOUPLTYPE(ir->etc),ETCOUPLTYPE(etcVRESCALE));
394 warning_note(wi,warn_buf);
397 if ((ir->etc==etcNOSEHOOVER || ir->etc==etcANDERSEN || ir->etc==etcANDERSENINTERVAL)
398 && ir->epc==epcBERENDSEN)
400 sprintf(warn_buf,"Using Berendsen pressure coupling invalidates the "
401 "true ensemble for the thermostat");
402 warning(wi,warn_buf);
405 /* PRESSURE COUPLING */
406 if (ir->epc == epcISOTROPIC)
408 ir->epc = epcBERENDSEN;
409 warning_note(wi,"Old option for pressure coupling given: "
410 "changing \"Isotropic\" to \"Berendsen\"\n");
413 if (ir->epc != epcNO)
415 dt_pcoupl = ir->nstpcouple*ir->delta_t;
417 sprintf(err_buf,"tau-p must be > 0 instead of %g\n",ir->tau_p);
418 CHECK(ir->tau_p <= 0);
420 if (ir->tau_p/dt_pcoupl < pcouple_min_integration_steps(ir->epc))
422 sprintf(warn_buf,"For proper integration of the %s barostat, tau-p (%g) should be at least %d times larger than nstpcouple*dt (%g)",
423 EPCOUPLTYPE(ir->epc),ir->tau_p,pcouple_min_integration_steps(ir->epc),dt_pcoupl);
424 warning(wi,warn_buf);
427 sprintf(err_buf,"compressibility must be > 0 when using pressure"
428 " coupling %s\n",EPCOUPLTYPE(ir->epc));
429 CHECK(ir->compress[XX][XX] < 0 || ir->compress[YY][YY] < 0 ||
430 ir->compress[ZZ][ZZ] < 0 ||
431 (trace(ir->compress) == 0 && ir->compress[YY][XX] <= 0 &&
432 ir->compress[ZZ][XX] <= 0 && ir->compress[ZZ][YY] <= 0));
434 sprintf(err_buf,"pressure coupling with PPPM not implemented, use PME");
435 CHECK(ir->coulombtype == eelPPPM);
438 else if (ir->coulombtype == eelPPPM)
440 sprintf(warn_buf,"The pressure with PPPM is incorrect, if you need the pressure use PME");
441 warning(wi,warn_buf);
448 if (ir->epc!=epcMTTK)
450 warning_error(wi,"NPT only defined for vv using Martyna-Tuckerman-Tobias-Klein equations");
456 /* More checks are in triple check (grompp.c) */
457 if (ir->coulombtype == eelPPPM)
459 warning_error(wi,"PPPM is not functional in the current version, we plan to implement PPPM through a small modification of the PME code");
462 if (ir->coulombtype == eelSWITCH) {
463 sprintf(warn_buf,"coulombtype = %s is only for testing purposes and can lead to serious artifacts, advice: use coulombtype = %s",
464 eel_names[ir->coulombtype],
465 eel_names[eelRF_ZERO]);
466 warning(wi,warn_buf);
469 if (ir->epsilon_r!=1 && ir->implicit_solvent==eisGBSA) {
470 sprintf(warn_buf,"epsilon-r = %g with GB implicit solvent, will use this value for inner dielectric",ir->epsilon_r);
471 warning_note(wi,warn_buf);
474 if (EEL_RF(ir->coulombtype) && ir->epsilon_rf==1 && ir->epsilon_r!=1) {
475 sprintf(warn_buf,"epsilon-r = %g and epsilon-rf = 1 with reaction field, assuming old format and exchanging epsilon-r and epsilon-rf",ir->epsilon_r);
476 warning(wi,warn_buf);
477 ir->epsilon_rf = ir->epsilon_r;
481 if (getenv("GALACTIC_DYNAMICS") == NULL) {
482 sprintf(err_buf,"epsilon-r must be >= 0 instead of %g\n",ir->epsilon_r);
483 CHECK(ir->epsilon_r < 0);
486 if (EEL_RF(ir->coulombtype)) {
487 /* reaction field (at the cut-off) */
489 if (ir->coulombtype == eelRF_ZERO) {
490 sprintf(err_buf,"With coulombtype = %s, epsilon-rf must be 0",
491 eel_names[ir->coulombtype]);
492 CHECK(ir->epsilon_rf != 0);
495 sprintf(err_buf,"epsilon-rf must be >= epsilon-r");
496 CHECK((ir->epsilon_rf < ir->epsilon_r && ir->epsilon_rf != 0) ||
497 (ir->epsilon_r == 0));
498 if (ir->epsilon_rf == ir->epsilon_r) {
499 sprintf(warn_buf,"Using epsilon-rf = epsilon-r with %s does not make sense",
500 eel_names[ir->coulombtype]);
501 warning(wi,warn_buf);
504 /* Allow rlist>rcoulomb for tabulated long range stuff. This just
505 * means the interaction is zero outside rcoulomb, but it helps to
506 * provide accurate energy conservation.
508 if (EEL_MIGHT_BE_ZERO_AT_CUTOFF(ir->coulombtype)) {
509 if (EEL_SWITCHED(ir->coulombtype)) {
511 "With coulombtype = %s rcoulomb_switch must be < rcoulomb",
512 eel_names[ir->coulombtype]);
513 CHECK(ir->rcoulomb_switch >= ir->rcoulomb);
515 } else if (ir->coulombtype == eelCUT || EEL_RF(ir->coulombtype)) {
516 sprintf(err_buf,"With coulombtype = %s, rcoulomb must be >= rlist",
517 eel_names[ir->coulombtype]);
518 CHECK(ir->rlist > ir->rcoulomb);
521 if (EEL_FULL(ir->coulombtype)) {
522 if (ir->coulombtype==eelPMESWITCH || ir->coulombtype==eelPMEUSER ||
523 ir->coulombtype==eelPMEUSERSWITCH) {
524 sprintf(err_buf,"With coulombtype = %s, rcoulomb must be <= rlist",
525 eel_names[ir->coulombtype]);
526 CHECK(ir->rcoulomb > ir->rlist);
528 if (ir->coulombtype == eelPME) {
530 "With coulombtype = %s, rcoulomb must be equal to rlist\n"
531 "If you want optimal energy conservation or exact integration use %s",
532 eel_names[ir->coulombtype],eel_names[eelPMESWITCH]);
535 "With coulombtype = %s, rcoulomb must be equal to rlist",
536 eel_names[ir->coulombtype]);
538 CHECK(ir->rcoulomb != ir->rlist);
542 if (EEL_PME(ir->coulombtype)) {
543 if (ir->pme_order < 3) {
544 warning_error(wi,"pme-order can not be smaller than 3");
548 if (ir->nwall==2 && EEL_FULL(ir->coulombtype)) {
549 if (ir->ewald_geometry == eewg3D) {
550 sprintf(warn_buf,"With pbc=%s you should use ewald-geometry=%s",
551 epbc_names[ir->ePBC],eewg_names[eewg3DC]);
552 warning(wi,warn_buf);
554 /* This check avoids extra pbc coding for exclusion corrections */
555 sprintf(err_buf,"wall-ewald-zfac should be >= 2");
556 CHECK(ir->wall_ewald_zfac < 2);
559 if (EVDW_SWITCHED(ir->vdwtype)) {
560 sprintf(err_buf,"With vdwtype = %s rvdw-switch must be < rvdw",
561 evdw_names[ir->vdwtype]);
562 CHECK(ir->rvdw_switch >= ir->rvdw);
563 } else if (ir->vdwtype == evdwCUT) {
564 sprintf(err_buf,"With vdwtype = %s, rvdw must be >= rlist",evdw_names[ir->vdwtype]);
565 CHECK(ir->rlist > ir->rvdw);
567 if (EEL_IS_ZERO_AT_CUTOFF(ir->coulombtype)
568 && (ir->rlistlong <= ir->rcoulomb)) {
569 sprintf(warn_buf,"For energy conservation with switch/shift potentials, %s should be 0.1 to 0.3 nm larger than rcoulomb.",
570 IR_TWINRANGE(*ir) ? "rlistlong" : "rlist");
571 warning_note(wi,warn_buf);
573 if (EVDW_SWITCHED(ir->vdwtype) && (ir->rlistlong <= ir->rvdw)) {
574 sprintf(warn_buf,"For energy conservation with switch/shift potentials, %s should be 0.1 to 0.3 nm larger than rvdw.",
575 IR_TWINRANGE(*ir) ? "rlistlong" : "rlist");
576 warning_note(wi,warn_buf);
579 if (ir->vdwtype == evdwUSER && ir->eDispCorr != edispcNO) {
580 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.");
583 if (ir->nstlist == -1) {
585 "nstlist=-1 only works with switched or shifted potentials,\n"
586 "suggestion: use vdw-type=%s and coulomb-type=%s",
587 evdw_names[evdwSHIFT],eel_names[eelPMESWITCH]);
588 CHECK(!(EEL_MIGHT_BE_ZERO_AT_CUTOFF(ir->coulombtype) &&
589 EVDW_MIGHT_BE_ZERO_AT_CUTOFF(ir->vdwtype)));
591 sprintf(err_buf,"With nstlist=-1 rvdw and rcoulomb should be smaller than rlist to account for diffusion and possibly charge-group radii");
592 CHECK(ir->rvdw >= ir->rlist || ir->rcoulomb >= ir->rlist);
594 sprintf(err_buf,"nstlist can not be smaller than -1");
595 CHECK(ir->nstlist < -1);
597 if (ir->eI == eiLBFGS && (ir->coulombtype==eelCUT || ir->vdwtype==evdwCUT)
599 warning(wi,"For efficient BFGS minimization, use switch/shift/pme instead of cut-off.");
602 if (ir->eI == eiLBFGS && ir->nbfgscorr <= 0) {
603 warning(wi,"Using L-BFGS with nbfgscorr<=0 just gets you steepest descent.");
607 if (ir->efep != efepNO) {
608 sprintf(err_buf,"The soft-core power is %d and can only be 1 or 2",
610 CHECK(ir->sc_alpha!=0 && ir->sc_power!=1 && ir->sc_power!=2);
613 /* ENERGY CONSERVATION */
616 if (!EVDW_MIGHT_BE_ZERO_AT_CUTOFF(ir->vdwtype) && ir->rvdw > 0)
618 sprintf(warn_buf,"You are using a cut-off for VdW interactions with NVE, for good energy conservation use vdwtype = %s (possibly with DispCorr)",
619 evdw_names[evdwSHIFT]);
620 warning_note(wi,warn_buf);
622 if (!EEL_MIGHT_BE_ZERO_AT_CUTOFF(ir->coulombtype) && ir->rcoulomb > 0)
624 sprintf(warn_buf,"You are using a cut-off for electrostatics with NVE, for good energy conservation use coulombtype = %s or %s",
625 eel_names[eelPMESWITCH],eel_names[eelRF_ZERO]);
626 warning_note(wi,warn_buf);
630 /* IMPLICIT SOLVENT */
631 if(ir->coulombtype==eelGB_NOTUSED)
633 ir->coulombtype=eelCUT;
634 ir->implicit_solvent=eisGBSA;
635 fprintf(stderr,"Note: Old option for generalized born electrostatics given:\n"
636 "Changing coulombtype from \"generalized-born\" to \"cut-off\" and instead\n"
637 "setting implicit-solvent value to \"GBSA\" in input section.\n");
640 if(ir->sa_algorithm==esaSTILL)
642 sprintf(err_buf,"Still SA algorithm not available yet, use %s or %s instead\n",esa_names[esaAPPROX],esa_names[esaNO]);
643 CHECK(ir->sa_algorithm == esaSTILL);
646 if(ir->implicit_solvent==eisGBSA)
648 sprintf(err_buf,"With GBSA implicit solvent, rgbradii must be equal to rlist.");
649 CHECK(ir->rgbradii != ir->rlist);
651 if(ir->coulombtype!=eelCUT)
653 sprintf(err_buf,"With GBSA, coulombtype must be equal to %s\n",eel_names[eelCUT]);
654 CHECK(ir->coulombtype!=eelCUT);
656 if(ir->vdwtype!=evdwCUT)
658 sprintf(err_buf,"With GBSA, vdw-type must be equal to %s\n",evdw_names[evdwCUT]);
659 CHECK(ir->vdwtype!=evdwCUT);
663 sprintf(warn_buf,"Using GBSA with nstgbradii<1, setting nstgbradii=1");
664 warning_note(wi,warn_buf);
667 if(ir->sa_algorithm==esaNO)
669 sprintf(warn_buf,"No SA (non-polar) calculation requested together with GB. Are you sure this is what you want?\n");
670 warning_note(wi,warn_buf);
672 if(ir->sa_surface_tension<0 && ir->sa_algorithm!=esaNO)
674 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");
675 warning_note(wi,warn_buf);
677 if(ir->gb_algorithm==egbSTILL)
679 ir->sa_surface_tension = 0.0049 * CAL2JOULE * 100;
683 ir->sa_surface_tension = 0.0054 * CAL2JOULE * 100;
686 if(ir->sa_surface_tension==0 && ir->sa_algorithm!=esaNO)
688 sprintf(err_buf, "Surface tension set to 0 while SA-calculation requested\n");
689 CHECK(ir->sa_surface_tension==0 && ir->sa_algorithm!=esaNO);
695 static int str_nelem(const char *str,int maxptr,char *ptr[])
703 while (*copy != '\0') {
705 gmx_fatal(FARGS,"Too many groups on line: '%s' (max is %d)",
710 while ((*copy != '\0') && !isspace(*copy))
724 static void parse_n_double(char *str,int *n,double **r)
729 *n = str_nelem(str,MAXPTR,ptr);
732 for(i=0; i<*n; i++) {
733 (*r)[i] = strtod(ptr[i],NULL);
737 static void do_wall_params(t_inputrec *ir,
738 char *wall_atomtype, char *wall_density,
745 opts->wall_atomtype[0] = NULL;
746 opts->wall_atomtype[1] = NULL;
748 ir->wall_atomtype[0] = -1;
749 ir->wall_atomtype[1] = -1;
750 ir->wall_density[0] = 0;
751 ir->wall_density[1] = 0;
755 nstr = str_nelem(wall_atomtype,MAXPTR,names);
756 if (nstr != ir->nwall)
758 gmx_fatal(FARGS,"Expected %d elements for wall_atomtype, found %d",
761 for(i=0; i<ir->nwall; i++)
763 opts->wall_atomtype[i] = strdup(names[i]);
766 if (ir->wall_type == ewt93 || ir->wall_type == ewt104) {
767 nstr = str_nelem(wall_density,MAXPTR,names);
768 if (nstr != ir->nwall)
770 gmx_fatal(FARGS,"Expected %d elements for wall-density, found %d",ir->nwall,nstr);
772 for(i=0; i<ir->nwall; i++)
774 sscanf(names[i],"%lf",&dbl);
777 gmx_fatal(FARGS,"wall-density[%d] = %f\n",i,dbl);
779 ir->wall_density[i] = dbl;
785 static void add_wall_energrps(gmx_groups_t *groups,int nwall,t_symtab *symtab)
792 srenew(groups->grpname,groups->ngrpname+nwall);
793 grps = &(groups->grps[egcENER]);
794 srenew(grps->nm_ind,grps->nr+nwall);
795 for(i=0; i<nwall; i++) {
796 sprintf(str,"wall%d",i);
797 groups->grpname[groups->ngrpname] = put_symtab(symtab,str);
798 grps->nm_ind[grps->nr++] = groups->ngrpname++;
803 void get_ir(const char *mdparin,const char *mdparout,
804 t_inputrec *ir,t_gromppopts *opts,
812 char warn_buf[STRLEN];
814 inp = read_inpfile(mdparin, &ninp, NULL, wi);
816 snew(dumstr[0],STRLEN);
817 snew(dumstr[1],STRLEN);
821 REM_TYPE("domain-decomposition");
822 REPL_TYPE("unconstrained-start","continuation");
823 REM_TYPE("dihre-tau");
824 REM_TYPE("nstdihreout");
825 REM_TYPE("nstcheckpoint");
827 CCTYPE ("VARIOUS PREPROCESSING OPTIONS");
828 CTYPE ("Preprocessor information: use cpp syntax.");
829 CTYPE ("e.g.: -I/home/joe/doe -I/home/mary/roe");
830 STYPE ("include", opts->include, NULL);
831 CTYPE ("e.g.: -DPOSRES -DFLEXIBLE (note these variable names are case sensitive)");
832 STYPE ("define", opts->define, NULL);
834 CCTYPE ("RUN CONTROL PARAMETERS");
835 EETYPE("integrator", ir->eI, ei_names);
836 CTYPE ("Start time and timestep in ps");
837 RTYPE ("tinit", ir->init_t, 0.0);
838 RTYPE ("dt", ir->delta_t, 0.001);
839 STEPTYPE ("nsteps", ir->nsteps, 0);
840 CTYPE ("For exact run continuation or redoing part of a run");
841 STEPTYPE ("init-step",ir->init_step, 0);
842 CTYPE ("Part index is updated automatically on checkpointing (keeps files separate)");
843 ITYPE ("simulation-part", ir->simulation_part, 1);
844 CTYPE ("mode for center of mass motion removal");
845 EETYPE("comm-mode", ir->comm_mode, ecm_names);
846 CTYPE ("number of steps for center of mass motion removal");
847 ITYPE ("nstcomm", ir->nstcomm, 10);
848 CTYPE ("group(s) for center of mass motion removal");
849 STYPE ("comm-grps", vcm, NULL);
851 CCTYPE ("LANGEVIN DYNAMICS OPTIONS");
852 CTYPE ("Friction coefficient (amu/ps) and random seed");
853 RTYPE ("bd-fric", ir->bd_fric, 0.0);
854 ITYPE ("ld-seed", ir->ld_seed, 1993);
857 CCTYPE ("ENERGY MINIMIZATION OPTIONS");
858 CTYPE ("Force tolerance and initial step-size");
859 RTYPE ("emtol", ir->em_tol, 10.0);
860 RTYPE ("emstep", ir->em_stepsize,0.01);
861 CTYPE ("Max number of iterations in relax-shells");
862 ITYPE ("niter", ir->niter, 20);
863 CTYPE ("Step size (ps^2) for minimization of flexible constraints");
864 RTYPE ("fcstep", ir->fc_stepsize, 0);
865 CTYPE ("Frequency of steepest descents steps when doing CG");
866 ITYPE ("nstcgsteep", ir->nstcgsteep, 1000);
867 ITYPE ("nbfgscorr", ir->nbfgscorr, 10);
869 CCTYPE ("TEST PARTICLE INSERTION OPTIONS");
870 RTYPE ("rtpi", ir->rtpi, 0.05);
873 CCTYPE ("OUTPUT CONTROL OPTIONS");
874 CTYPE ("Output frequency for coords (x), velocities (v) and forces (f)");
875 ITYPE ("nstxout", ir->nstxout, 100);
876 ITYPE ("nstvout", ir->nstvout, 100);
877 ITYPE ("nstfout", ir->nstfout, 0);
878 ir->nstcheckpoint = 1000;
879 CTYPE ("Output frequency for energies to log file and energy file");
880 ITYPE ("nstlog", ir->nstlog, 100);
881 ITYPE ("nstcalcenergy",ir->nstcalcenergy, -1);
882 ITYPE ("nstenergy", ir->nstenergy, 100);
883 CTYPE ("Output frequency and precision for .xtc file");
884 ITYPE ("nstxtcout", ir->nstxtcout, 0);
885 RTYPE ("xtc-precision",ir->xtcprec, 1000.0);
886 CTYPE ("This selects the subset of atoms for the .xtc file. You can");
887 CTYPE ("select multiple groups. By default all atoms will be written.");
888 STYPE ("xtc-grps", xtc_grps, NULL);
889 CTYPE ("Selection of energy groups");
890 STYPE ("energygrps", energy, NULL);
892 /* Neighbor searching */
893 CCTYPE ("NEIGHBORSEARCHING PARAMETERS");
894 CTYPE ("nblist update frequency");
895 ITYPE ("nstlist", ir->nstlist, 10);
896 CTYPE ("ns algorithm (simple or grid)");
897 EETYPE("ns-type", ir->ns_type, ens_names);
898 /* set ndelta to the optimal value of 2 */
900 CTYPE ("Periodic boundary conditions: xyz, no, xy");
901 EETYPE("pbc", ir->ePBC, epbc_names);
902 EETYPE("periodic-molecules", ir->bPeriodicMols, yesno_names);
903 CTYPE ("nblist cut-off");
904 RTYPE ("rlist", ir->rlist, 1.0);
905 CTYPE ("long-range cut-off for switched potentials");
906 RTYPE ("rlistlong", ir->rlistlong, -1);
909 CCTYPE ("OPTIONS FOR ELECTROSTATICS AND VDW");
910 CTYPE ("Method for doing electrostatics");
911 EETYPE("coulombtype", ir->coulombtype, eel_names);
912 CTYPE ("cut-off lengths");
913 RTYPE ("rcoulomb-switch", ir->rcoulomb_switch, 0.0);
914 RTYPE ("rcoulomb", ir->rcoulomb, 1.0);
915 CTYPE ("Relative dielectric constant for the medium and the reaction field");
916 RTYPE ("epsilon-r", ir->epsilon_r, 1.0);
917 RTYPE ("epsilon-rf", ir->epsilon_rf, 1.0);
918 CTYPE ("Method for doing Van der Waals");
919 EETYPE("vdw-type", ir->vdwtype, evdw_names);
920 CTYPE ("cut-off lengths");
921 RTYPE ("rvdw-switch", ir->rvdw_switch, 0.0);
922 RTYPE ("rvdw", ir->rvdw, 1.0);
923 CTYPE ("Apply long range dispersion corrections for Energy and Pressure");
924 EETYPE("DispCorr", ir->eDispCorr, edispc_names);
925 CTYPE ("Extension of the potential lookup tables beyond the cut-off");
926 RTYPE ("table-extension", ir->tabext, 1.0);
927 CTYPE ("Seperate tables between energy group pairs");
928 STYPE ("energygrp-table", egptable, NULL);
929 CTYPE ("Spacing for the PME/PPPM FFT grid");
930 RTYPE ("fourierspacing", opts->fourierspacing,0.12);
931 CTYPE ("FFT grid size, when a value is 0 fourierspacing will be used");
932 ITYPE ("fourier-nx", ir->nkx, 0);
933 ITYPE ("fourier-ny", ir->nky, 0);
934 ITYPE ("fourier-nz", ir->nkz, 0);
935 CTYPE ("EWALD/PME/PPPM parameters");
936 ITYPE ("pme-order", ir->pme_order, 4);
937 RTYPE ("ewald-rtol", ir->ewald_rtol, 0.00001);
938 EETYPE("ewald-geometry", ir->ewald_geometry, eewg_names);
939 RTYPE ("epsilon-surface", ir->epsilon_surface, 0.0);
940 EETYPE("optimize-fft",ir->bOptFFT, yesno_names);
942 CCTYPE("IMPLICIT SOLVENT ALGORITHM");
943 EETYPE("implicit-solvent", ir->implicit_solvent, eis_names);
945 CCTYPE ("GENERALIZED BORN ELECTROSTATICS");
946 CTYPE ("Algorithm for calculating Born radii");
947 EETYPE("gb-algorithm", ir->gb_algorithm, egb_names);
948 CTYPE ("Frequency of calculating the Born radii inside rlist");
949 ITYPE ("nstgbradii", ir->nstgbradii, 1);
950 CTYPE ("Cutoff for Born radii calculation; the contribution from atoms");
951 CTYPE ("between rlist and rgbradii is updated every nstlist steps");
952 RTYPE ("rgbradii", ir->rgbradii, 1.0);
953 CTYPE ("Dielectric coefficient of the implicit solvent");
954 RTYPE ("gb-epsilon-solvent",ir->gb_epsilon_solvent, 80.0);
955 CTYPE ("Salt concentration in M for Generalized Born models");
956 RTYPE ("gb-saltconc", ir->gb_saltconc, 0.0);
957 CTYPE ("Scaling factors used in the OBC GB model. Default values are OBC(II)");
958 RTYPE ("gb-obc-alpha", ir->gb_obc_alpha, 1.0);
959 RTYPE ("gb-obc-beta", ir->gb_obc_beta, 0.8);
960 RTYPE ("gb-obc-gamma", ir->gb_obc_gamma, 4.85);
961 RTYPE ("gb-dielectric-offset", ir->gb_dielectric_offset, 0.009);
962 EETYPE("sa-algorithm", ir->sa_algorithm, esa_names);
963 CTYPE ("Surface tension (kJ/mol/nm^2) for the SA (nonpolar surface) part of GBSA");
964 CTYPE ("The value -1 will set default value for Still/HCT/OBC GB-models.");
965 RTYPE ("sa-surface-tension", ir->sa_surface_tension, -1);
968 CCTYPE ("OPTIONS FOR WEAK COUPLING ALGORITHMS");
969 CTYPE ("Temperature coupling");
970 EETYPE("tcoupl", ir->etc, etcoupl_names);
971 ITYPE ("nsttcouple", ir->nsttcouple, -1);
972 ITYPE("nh-chain-length", ir->opts.nhchainlength, NHCHAINLENGTH);
973 CTYPE ("Groups to couple separately");
974 STYPE ("tc-grps", tcgrps, NULL);
975 CTYPE ("Time constant (ps) and reference temperature (K)");
976 STYPE ("tau-t", tau_t, NULL);
977 STYPE ("ref-t", ref_t, NULL);
978 CTYPE ("pressure coupling");
979 EETYPE("pcoupl", ir->epc, epcoupl_names);
980 EETYPE("pcoupltype", ir->epct, epcoupltype_names);
981 ITYPE ("nstpcouple", ir->nstpcouple, -1);
982 CTYPE ("Time constant (ps), compressibility (1/bar) and reference P (bar)");
983 RTYPE ("tau-p", ir->tau_p, 1.0);
984 STYPE ("compressibility", dumstr[0], NULL);
985 STYPE ("ref-p", dumstr[1], NULL);
986 CTYPE ("Scaling of reference coordinates, No, All or COM");
987 EETYPE ("refcoord-scaling",ir->refcoord_scaling,erefscaling_names);
989 CTYPE ("Random seed for Andersen thermostat");
990 ITYPE ("andersen-seed", ir->andersen_seed, 815131);
993 CCTYPE ("OPTIONS FOR QMMM calculations");
994 EETYPE("QMMM", ir->bQMMM, yesno_names);
995 CTYPE ("Groups treated Quantum Mechanically");
996 STYPE ("QMMM-grps", QMMM, NULL);
998 STYPE("QMmethod", QMmethod, NULL);
999 CTYPE ("QMMM scheme");
1000 EETYPE("QMMMscheme", ir->QMMMscheme, eQMMMscheme_names);
1001 CTYPE ("QM basisset");
1002 STYPE("QMbasis", QMbasis, NULL);
1003 CTYPE ("QM charge");
1004 STYPE ("QMcharge", QMcharge,NULL);
1005 CTYPE ("QM multiplicity");
1006 STYPE ("QMmult", QMmult,NULL);
1007 CTYPE ("Surface Hopping");
1008 STYPE ("SH", bSH, NULL);
1009 CTYPE ("CAS space options");
1010 STYPE ("CASorbitals", CASorbitals, NULL);
1011 STYPE ("CASelectrons", CASelectrons, NULL);
1012 STYPE ("SAon", SAon, NULL);
1013 STYPE ("SAoff",SAoff,NULL);
1014 STYPE ("SAsteps", SAsteps, NULL);
1015 CTYPE ("Scale factor for MM charges");
1016 RTYPE ("MMChargeScaleFactor", ir->scalefactor, 1.0);
1017 CTYPE ("Optimization of QM subsystem");
1018 STYPE ("bOPT", bOPT, NULL);
1019 STYPE ("bTS", bTS, NULL);
1021 /* Simulated annealing */
1022 CCTYPE("SIMULATED ANNEALING");
1023 CTYPE ("Type of annealing for each temperature group (no/single/periodic)");
1024 STYPE ("annealing", anneal, NULL);
1025 CTYPE ("Number of time points to use for specifying annealing in each group");
1026 STYPE ("annealing-npoints", anneal_npoints, NULL);
1027 CTYPE ("List of times at the annealing points for each group");
1028 STYPE ("annealing-time", anneal_time, NULL);
1029 CTYPE ("Temp. at each annealing point, for each group.");
1030 STYPE ("annealing-temp", anneal_temp, NULL);
1033 CCTYPE ("GENERATE VELOCITIES FOR STARTUP RUN");
1034 EETYPE("gen-vel", opts->bGenVel, yesno_names);
1035 RTYPE ("gen-temp", opts->tempi, 300.0);
1036 ITYPE ("gen-seed", opts->seed, 173529);
1039 CCTYPE ("OPTIONS FOR BONDS");
1040 EETYPE("constraints", opts->nshake, constraints);
1041 CTYPE ("Type of constraint algorithm");
1042 EETYPE("constraint-algorithm", ir->eConstrAlg, econstr_names);
1043 CTYPE ("Do not constrain the start configuration");
1044 EETYPE("continuation", ir->bContinuation, yesno_names);
1045 CTYPE ("Use successive overrelaxation to reduce the number of shake iterations");
1046 EETYPE("Shake-SOR", ir->bShakeSOR, yesno_names);
1047 CTYPE ("Relative tolerance of shake");
1048 RTYPE ("shake-tol", ir->shake_tol, 0.0001);
1049 CTYPE ("Highest order in the expansion of the constraint coupling matrix");
1050 ITYPE ("lincs-order", ir->nProjOrder, 4);
1051 CTYPE ("Number of iterations in the final step of LINCS. 1 is fine for");
1052 CTYPE ("normal simulations, but use 2 to conserve energy in NVE runs.");
1053 CTYPE ("For energy minimization with constraints it should be 4 to 8.");
1054 ITYPE ("lincs-iter", ir->nLincsIter, 1);
1055 CTYPE ("Lincs will write a warning to the stderr if in one step a bond");
1056 CTYPE ("rotates over more degrees than");
1057 RTYPE ("lincs-warnangle", ir->LincsWarnAngle, 30.0);
1058 CTYPE ("Convert harmonic bonds to morse potentials");
1059 EETYPE("morse", opts->bMorse,yesno_names);
1061 /* Energy group exclusions */
1062 CCTYPE ("ENERGY GROUP EXCLUSIONS");
1063 CTYPE ("Pairs of energy groups for which all non-bonded interactions are excluded");
1064 STYPE ("energygrp-excl", egpexcl, NULL);
1068 CTYPE ("Number of walls, type, atom types, densities and box-z scale factor for Ewald");
1069 ITYPE ("nwall", ir->nwall, 0);
1070 EETYPE("wall-type", ir->wall_type, ewt_names);
1071 RTYPE ("wall-r-linpot", ir->wall_r_linpot, -1);
1072 STYPE ("wall-atomtype", wall_atomtype, NULL);
1073 STYPE ("wall-density", wall_density, NULL);
1074 RTYPE ("wall-ewald-zfac", ir->wall_ewald_zfac, 3);
1077 CCTYPE("COM PULLING");
1078 CTYPE("Pull type: no, umbrella, constraint or constant-force");
1079 EETYPE("pull", ir->ePull, epull_names);
1080 if (ir->ePull != epullNO) {
1082 pull_grp = read_pullparams(&ninp,&inp,ir->pull,&opts->pull_start,wi);
1086 CCTYPE("NMR refinement stuff");
1087 CTYPE ("Distance restraints type: No, Simple or Ensemble");
1088 EETYPE("disre", ir->eDisre, edisre_names);
1089 CTYPE ("Force weighting of pairs in one distance restraint: Conservative or Equal");
1090 EETYPE("disre-weighting", ir->eDisreWeighting, edisreweighting_names);
1091 CTYPE ("Use sqrt of the time averaged times the instantaneous violation");
1092 EETYPE("disre-mixed", ir->bDisreMixed, yesno_names);
1093 RTYPE ("disre-fc", ir->dr_fc, 1000.0);
1094 RTYPE ("disre-tau", ir->dr_tau, 0.0);
1095 CTYPE ("Output frequency for pair distances to energy file");
1096 ITYPE ("nstdisreout", ir->nstdisreout, 100);
1097 CTYPE ("Orientation restraints: No or Yes");
1098 EETYPE("orire", opts->bOrire, yesno_names);
1099 CTYPE ("Orientation restraints force constant and tau for time averaging");
1100 RTYPE ("orire-fc", ir->orires_fc, 0.0);
1101 RTYPE ("orire-tau", ir->orires_tau, 0.0);
1102 STYPE ("orire-fitgrp",orirefitgrp, NULL);
1103 CTYPE ("Output frequency for trace(SD) and S to energy file");
1104 ITYPE ("nstorireout", ir->nstorireout, 100);
1105 CTYPE ("Dihedral angle restraints: No or Yes");
1106 EETYPE("dihre", opts->bDihre, yesno_names);
1107 RTYPE ("dihre-fc", ir->dihre_fc, 1000.0);
1109 /* Free energy stuff */
1110 CCTYPE ("Free energy control stuff");
1111 EETYPE("free-energy", ir->efep, efep_names);
1112 RTYPE ("init-lambda", ir->init_lambda,0.0);
1113 RTYPE ("delta-lambda",ir->delta_lambda,0.0);
1114 STYPE ("foreign-lambda", foreign_lambda, NULL);
1115 RTYPE ("sc-alpha",ir->sc_alpha,0.0);
1116 ITYPE ("sc-power",ir->sc_power,0);
1117 RTYPE ("sc-sigma",ir->sc_sigma,0.3);
1118 ITYPE ("nstdhdl", ir->nstdhdl, 10);
1119 EETYPE("separate-dhdl-file", ir->separate_dhdl_file,
1120 separate_dhdl_file_names);
1121 EETYPE("dhdl-derivatives", ir->dhdl_derivatives, dhdl_derivatives_names);
1122 ITYPE ("dh-hist-size", ir->dh_hist_size, 0);
1123 RTYPE ("dh-hist-spacing", ir->dh_hist_spacing, 0.1);
1124 STYPE ("couple-moltype", couple_moltype, NULL);
1125 EETYPE("couple-lambda0", opts->couple_lam0, couple_lam);
1126 EETYPE("couple-lambda1", opts->couple_lam1, couple_lam);
1127 EETYPE("couple-intramol", opts->bCoupleIntra, yesno_names);
1129 /* Non-equilibrium MD stuff */
1130 CCTYPE("Non-equilibrium MD stuff");
1131 STYPE ("acc-grps", accgrps, NULL);
1132 STYPE ("accelerate", acc, NULL);
1133 STYPE ("freezegrps", freeze, NULL);
1134 STYPE ("freezedim", frdim, NULL);
1135 RTYPE ("cos-acceleration", ir->cos_accel, 0);
1136 STYPE ("deform", deform, NULL);
1138 /* Electric fields */
1139 CCTYPE("Electric fields");
1140 CTYPE ("Format is number of terms (int) and for all terms an amplitude (real)");
1141 CTYPE ("and a phase angle (real)");
1142 STYPE ("E-x", efield_x, NULL);
1143 STYPE ("E-xt", efield_xt, NULL);
1144 STYPE ("E-y", efield_y, NULL);
1145 STYPE ("E-yt", efield_yt, NULL);
1146 STYPE ("E-z", efield_z, NULL);
1147 STYPE ("E-zt", efield_zt, NULL);
1149 /* User defined thingies */
1150 CCTYPE ("User defined thingies");
1151 STYPE ("user1-grps", user1, NULL);
1152 STYPE ("user2-grps", user2, NULL);
1153 ITYPE ("userint1", ir->userint1, 0);
1154 ITYPE ("userint2", ir->userint2, 0);
1155 ITYPE ("userint3", ir->userint3, 0);
1156 ITYPE ("userint4", ir->userint4, 0);
1157 RTYPE ("userreal1", ir->userreal1, 0);
1158 RTYPE ("userreal2", ir->userreal2, 0);
1159 RTYPE ("userreal3", ir->userreal3, 0);
1160 RTYPE ("userreal4", ir->userreal4, 0);
1163 write_inpfile(mdparout,ninp,inp,FALSE,wi);
1164 for (i=0; (i<ninp); i++) {
1166 sfree(inp[i].value);
1170 /* Process options if necessary */
1171 for(m=0; m<2; m++) {
1172 for(i=0; i<2*DIM; i++)
1177 if (sscanf(dumstr[m],"%lf",&(dumdub[m][XX]))!=1) {
1178 warning_error(wi,"Pressure coupling not enough values (I need 1)");
1180 dumdub[m][YY]=dumdub[m][ZZ]=dumdub[m][XX];
1182 case epctSEMIISOTROPIC:
1183 case epctSURFACETENSION:
1184 if (sscanf(dumstr[m],"%lf%lf",
1185 &(dumdub[m][XX]),&(dumdub[m][ZZ]))!=2) {
1186 warning_error(wi,"Pressure coupling not enough values (I need 2)");
1188 dumdub[m][YY]=dumdub[m][XX];
1190 case epctANISOTROPIC:
1191 if (sscanf(dumstr[m],"%lf%lf%lf%lf%lf%lf",
1192 &(dumdub[m][XX]),&(dumdub[m][YY]),&(dumdub[m][ZZ]),
1193 &(dumdub[m][3]),&(dumdub[m][4]),&(dumdub[m][5]))!=6) {
1194 warning_error(wi,"Pressure coupling not enough values (I need 6)");
1198 gmx_fatal(FARGS,"Pressure coupling type %s not implemented yet",
1199 epcoupltype_names[ir->epct]);
1203 clear_mat(ir->ref_p);
1204 clear_mat(ir->compress);
1205 for(i=0; i<DIM; i++) {
1206 ir->ref_p[i][i] = dumdub[1][i];
1207 ir->compress[i][i] = dumdub[0][i];
1209 if (ir->epct == epctANISOTROPIC) {
1210 ir->ref_p[XX][YY] = dumdub[1][3];
1211 ir->ref_p[XX][ZZ] = dumdub[1][4];
1212 ir->ref_p[YY][ZZ] = dumdub[1][5];
1213 if (ir->ref_p[XX][YY]!=0 && ir->ref_p[XX][ZZ]!=0 && ir->ref_p[YY][ZZ]!=0) {
1214 warning(wi,"All off-diagonal reference pressures are non-zero. Are you sure you want to apply a threefold shear stress?\n");
1216 ir->compress[XX][YY] = dumdub[0][3];
1217 ir->compress[XX][ZZ] = dumdub[0][4];
1218 ir->compress[YY][ZZ] = dumdub[0][5];
1219 for(i=0; i<DIM; i++) {
1220 for(m=0; m<i; m++) {
1221 ir->ref_p[i][m] = ir->ref_p[m][i];
1222 ir->compress[i][m] = ir->compress[m][i];
1227 if (ir->comm_mode == ecmNO)
1230 opts->couple_moltype = NULL;
1231 if (strlen(couple_moltype) > 0) {
1232 if (ir->efep != efepNO) {
1233 opts->couple_moltype = strdup(couple_moltype);
1234 if (opts->couple_lam0 == opts->couple_lam1)
1235 warning(wi,"The lambda=0 and lambda=1 states for coupling are identical");
1236 if (ir->eI == eiMD && (opts->couple_lam0 == ecouplamNONE ||
1237 opts->couple_lam1 == ecouplamNONE)) {
1238 warning(wi,"For proper sampling of the (nearly) decoupled state, stochastic dynamics should be used");
1241 warning(wi,"Can not couple a molecule with free-energy = no");
1245 do_wall_params(ir,wall_atomtype,wall_density,opts);
1247 if (opts->bOrire && str_nelem(orirefitgrp,MAXPTR,NULL)!=1) {
1248 warning_error(wi,"ERROR: Need one orientation restraint fit group\n");
1251 clear_mat(ir->deform);
1254 m = sscanf(deform,"%lf %lf %lf %lf %lf %lf",
1255 &(dumdub[0][0]),&(dumdub[0][1]),&(dumdub[0][2]),
1256 &(dumdub[0][3]),&(dumdub[0][4]),&(dumdub[0][5]));
1258 ir->deform[i][i] = dumdub[0][i];
1259 ir->deform[YY][XX] = dumdub[0][3];
1260 ir->deform[ZZ][XX] = dumdub[0][4];
1261 ir->deform[ZZ][YY] = dumdub[0][5];
1262 if (ir->epc != epcNO) {
1265 if (ir->deform[i][j]!=0 && ir->compress[i][j]!=0) {
1266 warning_error(wi,"A box element has deform set and compressibility > 0");
1270 if (ir->deform[i][j]!=0) {
1271 for(m=j; m<DIM; m++)
1272 if (ir->compress[m][j]!=0) {
1273 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.");
1274 warning(wi,warn_buf);
1279 if (ir->efep != efepNO) {
1280 parse_n_double(foreign_lambda,&ir->n_flambda,&ir->flambda);
1281 if (ir->n_flambda > 0 && ir->rlist < max(ir->rvdw,ir->rcoulomb)) {
1282 warning_note(wi,"For foreign lambda free energy differences it is assumed that the soft-core interactions have no effect beyond the neighborlist cut-off");
1292 static int search_QMstring(char *s,int ng,const char *gn[])
1294 /* same as normal search_string, but this one searches QM strings */
1297 for(i=0; (i<ng); i++)
1298 if (gmx_strcasecmp(s,gn[i]) == 0)
1301 gmx_fatal(FARGS,"this QM method or basisset (%s) is not implemented\n!",s);
1305 } /* search_QMstring */
1308 int search_string(char *s,int ng,char *gn[])
1312 for(i=0; (i<ng); i++)
1314 if (gmx_strcasecmp(s,gn[i]) == 0)
1320 gmx_fatal(FARGS,"Group %s not found in index file.\nGroup names must match either [moleculetype] names\nor custom index group names,in which case you\nmust supply an index file to the '-n' option of grompp.",s);
1325 static gmx_bool do_numbering(int natoms,gmx_groups_t *groups,int ng,char *ptrs[],
1326 t_blocka *block,char *gnames[],
1327 int gtype,int restnm,
1328 int grptp,gmx_bool bVerbose,
1331 unsigned short *cbuf;
1332 t_grps *grps=&(groups->grps[gtype]);
1333 int i,j,gid,aj,ognr,ntot=0;
1336 char warn_buf[STRLEN];
1340 fprintf(debug,"Starting numbering %d groups of type %d\n",ng,gtype);
1343 title = gtypes[gtype];
1346 /* Mark all id's as not set */
1347 for(i=0; (i<natoms); i++)
1352 snew(grps->nm_ind,ng+1); /* +1 for possible rest group */
1353 for(i=0; (i<ng); i++)
1355 /* Lookup the group name in the block structure */
1356 gid = search_string(ptrs[i],block->nr,gnames);
1357 if ((grptp != egrptpONE) || (i == 0))
1359 grps->nm_ind[grps->nr++]=gid;
1363 fprintf(debug,"Found gid %d for group %s\n",gid,ptrs[i]);
1366 /* Now go over the atoms in the group */
1367 for(j=block->index[gid]; (j<block->index[gid+1]); j++)
1372 /* Range checking */
1373 if ((aj < 0) || (aj >= natoms))
1375 gmx_fatal(FARGS,"Invalid atom number %d in indexfile",aj);
1377 /* Lookup up the old group number */
1381 gmx_fatal(FARGS,"Atom %d in multiple %s groups (%d and %d)",
1382 aj+1,title,ognr+1,i+1);
1386 /* Store the group number in buffer */
1387 if (grptp == egrptpONE)
1400 /* Now check whether we have done all atoms */
1404 if (grptp == egrptpALL)
1406 gmx_fatal(FARGS,"%d atoms are not part of any of the %s groups",
1409 else if (grptp == egrptpPART)
1411 sprintf(warn_buf,"%d atoms are not part of any of the %s groups",
1413 warning_note(wi,warn_buf);
1415 /* Assign all atoms currently unassigned to a rest group */
1416 for(j=0; (j<natoms); j++)
1418 if (cbuf[j] == NOGID)
1424 if (grptp != egrptpPART)
1429 "Making dummy/rest group for %s containing %d elements\n",
1432 /* Add group name "rest" */
1433 grps->nm_ind[grps->nr] = restnm;
1435 /* Assign the rest name to all atoms not currently assigned to a group */
1436 for(j=0; (j<natoms); j++)
1438 if (cbuf[j] == NOGID)
1449 groups->ngrpnr[gtype] = 0;
1450 groups->grpnr[gtype] = NULL;
1454 groups->ngrpnr[gtype] = natoms;
1455 snew(groups->grpnr[gtype],natoms);
1456 for(j=0; (j<natoms); j++)
1458 groups->grpnr[gtype][j] = cbuf[j];
1464 return (bRest && grptp == egrptpPART);
1467 static void calc_nrdf(gmx_mtop_t *mtop,t_inputrec *ir,char **gnames)
1470 gmx_groups_t *groups;
1472 int natoms,ai,aj,i,j,d,g,imin,jmin,nc;
1474 int *nrdf2,*na_vcm,na_tot;
1475 double *nrdf_tc,*nrdf_vcm,nrdf_uc,n_sub=0;
1476 gmx_mtop_atomloop_all_t aloop;
1478 int mb,mol,ftype,as;
1479 gmx_molblock_t *molb;
1480 gmx_moltype_t *molt;
1483 * First calc 3xnr-atoms for each group
1484 * then subtract half a degree of freedom for each constraint
1486 * Only atoms and nuclei contribute to the degrees of freedom...
1491 groups = &mtop->groups;
1492 natoms = mtop->natoms;
1494 /* Allocate one more for a possible rest group */
1495 /* We need to sum degrees of freedom into doubles,
1496 * since floats give too low nrdf's above 3 million atoms.
1498 snew(nrdf_tc,groups->grps[egcTC].nr+1);
1499 snew(nrdf_vcm,groups->grps[egcVCM].nr+1);
1500 snew(na_vcm,groups->grps[egcVCM].nr+1);
1502 for(i=0; i<groups->grps[egcTC].nr; i++)
1504 for(i=0; i<groups->grps[egcVCM].nr+1; i++)
1508 aloop = gmx_mtop_atomloop_all_init(mtop);
1509 while (gmx_mtop_atomloop_all_next(aloop,&i,&atom)) {
1511 if (atom->ptype == eptAtom || atom->ptype == eptNucleus) {
1512 g = ggrpnr(groups,egcFREEZE,i);
1513 /* Double count nrdf for particle i */
1514 for(d=0; d<DIM; d++) {
1515 if (opts->nFreeze[g][d] == 0) {
1519 nrdf_tc [ggrpnr(groups,egcTC ,i)] += 0.5*nrdf2[i];
1520 nrdf_vcm[ggrpnr(groups,egcVCM,i)] += 0.5*nrdf2[i];
1525 for(mb=0; mb<mtop->nmolblock; mb++) {
1526 molb = &mtop->molblock[mb];
1527 molt = &mtop->moltype[molb->type];
1528 atom = molt->atoms.atom;
1529 for(mol=0; mol<molb->nmol; mol++) {
1530 for (ftype=F_CONSTR; ftype<=F_CONSTRNC; ftype++) {
1531 ia = molt->ilist[ftype].iatoms;
1532 for(i=0; i<molt->ilist[ftype].nr; ) {
1533 /* Subtract degrees of freedom for the constraints,
1534 * if the particles still have degrees of freedom left.
1535 * If one of the particles is a vsite or a shell, then all
1536 * constraint motion will go there, but since they do not
1537 * contribute to the constraints the degrees of freedom do not
1542 if (((atom[ia[1]].ptype == eptNucleus) ||
1543 (atom[ia[1]].ptype == eptAtom)) &&
1544 ((atom[ia[2]].ptype == eptNucleus) ||
1545 (atom[ia[2]].ptype == eptAtom))) {
1554 imin = min(imin,nrdf2[ai]);
1555 jmin = min(jmin,nrdf2[aj]);
1558 nrdf_tc [ggrpnr(groups,egcTC ,ai)] -= 0.5*imin;
1559 nrdf_tc [ggrpnr(groups,egcTC ,aj)] -= 0.5*jmin;
1560 nrdf_vcm[ggrpnr(groups,egcVCM,ai)] -= 0.5*imin;
1561 nrdf_vcm[ggrpnr(groups,egcVCM,aj)] -= 0.5*jmin;
1563 ia += interaction_function[ftype].nratoms+1;
1564 i += interaction_function[ftype].nratoms+1;
1567 ia = molt->ilist[F_SETTLE].iatoms;
1568 for(i=0; i<molt->ilist[F_SETTLE].nr; ) {
1569 /* Subtract 1 dof from every atom in the SETTLE */
1570 for(ai=as+ia[1]; ai<as+ia[1]+3; ai++) {
1571 imin = min(2,nrdf2[ai]);
1573 nrdf_tc [ggrpnr(groups,egcTC ,ai)] -= 0.5*imin;
1574 nrdf_vcm[ggrpnr(groups,egcVCM,ai)] -= 0.5*imin;
1579 as += molt->atoms.nr;
1583 if (ir->ePull == epullCONSTRAINT) {
1584 /* Correct nrdf for the COM constraints.
1585 * We correct using the TC and VCM group of the first atom
1586 * in the reference and pull group. If atoms in one pull group
1587 * belong to different TC or VCM groups it is anyhow difficult
1588 * to determine the optimal nrdf assignment.
1591 if (pull->eGeom == epullgPOS) {
1593 for(i=0; i<DIM; i++) {
1600 for(i=0; i<pull->ngrp; i++) {
1602 if (pull->grp[0].nat > 0) {
1603 /* Subtract 1/2 dof from the reference group */
1604 ai = pull->grp[0].ind[0];
1605 if (nrdf_tc[ggrpnr(groups,egcTC,ai)] > 1) {
1606 nrdf_tc [ggrpnr(groups,egcTC ,ai)] -= 0.5;
1607 nrdf_vcm[ggrpnr(groups,egcVCM,ai)] -= 0.5;
1611 /* Subtract 1/2 dof from the pulled group */
1612 ai = pull->grp[1+i].ind[0];
1613 nrdf_tc [ggrpnr(groups,egcTC ,ai)] -= 0.5*imin;
1614 nrdf_vcm[ggrpnr(groups,egcVCM,ai)] -= 0.5*imin;
1615 if (nrdf_tc[ggrpnr(groups,egcTC,ai)] < 0)
1616 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)]]);
1620 if (ir->nstcomm != 0) {
1621 /* Subtract 3 from the number of degrees of freedom in each vcm group
1622 * when com translation is removed and 6 when rotation is removed
1625 switch (ir->comm_mode) {
1627 n_sub = ndof_com(ir);
1634 gmx_incons("Checking comm_mode");
1637 for(i=0; i<groups->grps[egcTC].nr; i++) {
1638 /* Count the number of atoms of TC group i for every VCM group */
1639 for(j=0; j<groups->grps[egcVCM].nr+1; j++)
1642 for(ai=0; ai<natoms; ai++)
1643 if (ggrpnr(groups,egcTC,ai) == i) {
1644 na_vcm[ggrpnr(groups,egcVCM,ai)]++;
1647 /* Correct for VCM removal according to the fraction of each VCM
1648 * group present in this TC group.
1650 nrdf_uc = nrdf_tc[i];
1652 fprintf(debug,"T-group[%d] nrdf_uc = %g, n_sub = %g\n",
1656 for(j=0; j<groups->grps[egcVCM].nr+1; j++) {
1657 if (nrdf_vcm[j] > n_sub) {
1658 nrdf_tc[i] += nrdf_uc*((double)na_vcm[j]/(double)na_tot)*
1659 (nrdf_vcm[j] - n_sub)/nrdf_vcm[j];
1662 fprintf(debug," nrdf_vcm[%d] = %g, nrdf = %g\n",
1663 j,nrdf_vcm[j],nrdf_tc[i]);
1668 for(i=0; (i<groups->grps[egcTC].nr); i++) {
1669 opts->nrdf[i] = nrdf_tc[i];
1670 if (opts->nrdf[i] < 0)
1673 "Number of degrees of freedom in T-Coupling group %s is %.2f\n",
1674 gnames[groups->grps[egcTC].nm_ind[i]],opts->nrdf[i]);
1683 static void decode_cos(char *s,t_cosines *cosine,gmx_bool bTime)
1686 char format[STRLEN],f1[STRLEN];
1697 sscanf(t,"%d",&(cosine->n));
1698 if (cosine->n <= 0) {
1701 snew(cosine->a,cosine->n);
1702 snew(cosine->phi,cosine->n);
1704 sprintf(format,"%%*d");
1705 for(i=0; (i<cosine->n); i++) {
1707 strcat(f1,"%lf%lf");
1708 if (sscanf(t,f1,&a,&phi) < 2)
1709 gmx_fatal(FARGS,"Invalid input for electric field shift: '%s'",t);
1712 strcat(format,"%*lf%*lf");
1719 static gmx_bool do_egp_flag(t_inputrec *ir,gmx_groups_t *groups,
1720 const char *option,const char *val,int flag)
1722 /* The maximum number of energy group pairs would be MAXPTR*(MAXPTR+1)/2.
1723 * But since this is much larger than STRLEN, such a line can not be parsed.
1724 * The real maximum is the number of names that fit in a string: STRLEN/2.
1726 #define EGP_MAX (STRLEN/2)
1728 char *names[EGP_MAX];
1732 gnames = groups->grpname;
1734 nelem = str_nelem(val,EGP_MAX,names);
1736 gmx_fatal(FARGS,"The number of groups for %s is odd",option);
1737 nr = groups->grps[egcENER].nr;
1739 for(i=0; i<nelem/2; i++) {
1742 gmx_strcasecmp(names[2*i],*(gnames[groups->grps[egcENER].nm_ind[j]])))
1745 gmx_fatal(FARGS,"%s in %s is not an energy group\n",
1749 gmx_strcasecmp(names[2*i+1],*(gnames[groups->grps[egcENER].nm_ind[k]])))
1752 gmx_fatal(FARGS,"%s in %s is not an energy group\n",
1753 names[2*i+1],option);
1754 if ((j < nr) && (k < nr)) {
1755 ir->opts.egp_flags[nr*j+k] |= flag;
1756 ir->opts.egp_flags[nr*k+j] |= flag;
1764 void do_index(const char* mdparin, const char *ndx,
1767 t_inputrec *ir,rvec *v,
1771 gmx_groups_t *groups;
1775 char warnbuf[STRLEN],**gnames;
1776 int nr,ntcg,ntau_t,nref_t,nacc,nofg,nSA,nSA_points,nSA_time,nSA_temp;
1779 int nacg,nfreeze,nfrdim,nenergy,nvcm,nuser;
1780 char *ptr1[MAXPTR],*ptr2[MAXPTR],*ptr3[MAXPTR];
1783 gmx_bool bExcl,bTable,bSetTCpar,bAnneal,bRest;
1784 int nQMmethod,nQMbasis,nQMcharge,nQMmult,nbSH,nCASorb,nCASelec,
1785 nSAon,nSAoff,nSAsteps,nQMg,nbOPT,nbTS;
1786 char warn_buf[STRLEN];
1789 fprintf(stderr,"processing index file...\n");
1793 snew(grps->index,1);
1795 atoms_all = gmx_mtop_global_atoms(mtop);
1796 analyse(&atoms_all,grps,&gnames,FALSE,TRUE);
1797 free_t_atoms(&atoms_all,FALSE);
1799 grps = init_index(ndx,&gnames);
1802 groups = &mtop->groups;
1803 natoms = mtop->natoms;
1804 symtab = &mtop->symtab;
1806 snew(groups->grpname,grps->nr+1);
1808 for(i=0; (i<grps->nr); i++) {
1809 groups->grpname[i] = put_symtab(symtab,gnames[i]);
1811 groups->grpname[i] = put_symtab(symtab,"rest");
1813 srenew(gnames,grps->nr+1);
1814 gnames[restnm] = *(groups->grpname[i]);
1815 groups->ngrpname = grps->nr+1;
1817 set_warning_line(wi,mdparin,-1);
1819 ntau_t = str_nelem(tau_t,MAXPTR,ptr1);
1820 nref_t = str_nelem(ref_t,MAXPTR,ptr2);
1821 ntcg = str_nelem(tcgrps,MAXPTR,ptr3);
1822 if ((ntau_t != ntcg) || (nref_t != ntcg)) {
1823 gmx_fatal(FARGS,"Invalid T coupling input: %d groups, %d ref-t values and "
1824 "%d tau-t values",ntcg,nref_t,ntau_t);
1827 bSetTCpar = (ir->etc || EI_SD(ir->eI) || ir->eI==eiBD || EI_TPI(ir->eI));
1828 do_numbering(natoms,groups,ntcg,ptr3,grps,gnames,egcTC,
1829 restnm,bSetTCpar ? egrptpALL : egrptpALL_GENREST,bVerbose,wi);
1830 nr = groups->grps[egcTC].nr;
1832 snew(ir->opts.nrdf,nr);
1833 snew(ir->opts.tau_t,nr);
1834 snew(ir->opts.ref_t,nr);
1835 if (ir->eI==eiBD && ir->bd_fric==0) {
1836 fprintf(stderr,"bd-fric=0, so tau-t will be used as the inverse friction constant(s)\n");
1843 gmx_fatal(FARGS,"Not enough ref-t and tau-t values!");
1847 for(i=0; (i<nr); i++)
1849 ir->opts.tau_t[i] = strtod(ptr1[i],NULL);
1850 if ((ir->eI == eiBD || ir->eI == eiSD2) && ir->opts.tau_t[i] <= 0)
1852 sprintf(warn_buf,"With integrator %s tau-t should be larger than 0",ei_names[ir->eI]);
1853 warning_error(wi,warn_buf);
1855 if ((ir->etc == etcVRESCALE && ir->opts.tau_t[i] >= 0) ||
1856 (ir->etc != etcVRESCALE && ir->opts.tau_t[i] > 0))
1858 tau_min = min(tau_min,ir->opts.tau_t[i]);
1861 if (ir->etc != etcNO && ir->nsttcouple == -1)
1863 ir->nsttcouple = ir_optimal_nsttcouple(ir);
1867 if ((ir->epc==epcMTTK) && (ir->etc>etcNO))
1870 mincouple = ir->nsttcouple;
1871 if (ir->nstpcouple < mincouple)
1873 mincouple = ir->nstpcouple;
1875 ir->nstpcouple = mincouple;
1876 ir->nsttcouple = mincouple;
1877 warning_note(wi,"for current Trotter decomposition methods with vv, nsttcouple and nstpcouple must be equal. Both have been reset to min(nsttcouple,nstpcouple)");
1880 nstcmin = tcouple_min_integration_steps(ir->etc);
1883 if (tau_min/(ir->delta_t*ir->nsttcouple) < nstcmin)
1885 sprintf(warn_buf,"For proper integration of the %s thermostat, tau-t (%g) should be at least %d times larger than nsttcouple*dt (%g)",
1886 ETCOUPLTYPE(ir->etc),
1888 ir->nsttcouple*ir->delta_t);
1889 warning(wi,warn_buf);
1892 for(i=0; (i<nr); i++)
1894 ir->opts.ref_t[i] = strtod(ptr2[i],NULL);
1895 if (ir->opts.ref_t[i] < 0)
1897 gmx_fatal(FARGS,"ref-t for group %d negative",i);
1902 /* Simulated annealing for each group. There are nr groups */
1903 nSA = str_nelem(anneal,MAXPTR,ptr1);
1904 if (nSA == 1 && (ptr1[0][0]=='n' || ptr1[0][0]=='N'))
1906 if(nSA>0 && nSA != nr)
1907 gmx_fatal(FARGS,"Not enough annealing values: %d (for %d groups)\n",nSA,nr);
1909 snew(ir->opts.annealing,nr);
1910 snew(ir->opts.anneal_npoints,nr);
1911 snew(ir->opts.anneal_time,nr);
1912 snew(ir->opts.anneal_temp,nr);
1914 ir->opts.annealing[i]=eannNO;
1915 ir->opts.anneal_npoints[i]=0;
1916 ir->opts.anneal_time[i]=NULL;
1917 ir->opts.anneal_temp[i]=NULL;
1922 if(ptr1[i][0]=='n' || ptr1[i][0]=='N') {
1923 ir->opts.annealing[i]=eannNO;
1924 } else if(ptr1[i][0]=='s'|| ptr1[i][0]=='S') {
1925 ir->opts.annealing[i]=eannSINGLE;
1927 } else if(ptr1[i][0]=='p'|| ptr1[i][0]=='P') {
1928 ir->opts.annealing[i]=eannPERIODIC;
1933 /* Read the other fields too */
1934 nSA_points = str_nelem(anneal_npoints,MAXPTR,ptr1);
1936 gmx_fatal(FARGS,"Found %d annealing-npoints values for %d groups\n",nSA_points,nSA);
1937 for(k=0,i=0;i<nr;i++) {
1938 ir->opts.anneal_npoints[i]=strtol(ptr1[i],NULL,10);
1939 if(ir->opts.anneal_npoints[i]==1)
1940 gmx_fatal(FARGS,"Please specify at least a start and an end point for annealing\n");
1941 snew(ir->opts.anneal_time[i],ir->opts.anneal_npoints[i]);
1942 snew(ir->opts.anneal_temp[i],ir->opts.anneal_npoints[i]);
1943 k += ir->opts.anneal_npoints[i];
1946 nSA_time = str_nelem(anneal_time,MAXPTR,ptr1);
1948 gmx_fatal(FARGS,"Found %d annealing-time values, wanter %d\n",nSA_time,k);
1949 nSA_temp = str_nelem(anneal_temp,MAXPTR,ptr2);
1951 gmx_fatal(FARGS,"Found %d annealing-temp values, wanted %d\n",nSA_temp,k);
1953 for(i=0,k=0;i<nr;i++) {
1955 for(j=0;j<ir->opts.anneal_npoints[i];j++) {
1956 ir->opts.anneal_time[i][j]=strtod(ptr1[k],NULL);
1957 ir->opts.anneal_temp[i][j]=strtod(ptr2[k],NULL);
1959 if(ir->opts.anneal_time[i][0] > (ir->init_t+GMX_REAL_EPS))
1960 gmx_fatal(FARGS,"First time point for annealing > init_t.\n");
1963 if(ir->opts.anneal_time[i][j]<ir->opts.anneal_time[i][j-1])
1964 gmx_fatal(FARGS,"Annealing timepoints out of order: t=%f comes after t=%f\n",
1965 ir->opts.anneal_time[i][j],ir->opts.anneal_time[i][j-1]);
1967 if(ir->opts.anneal_temp[i][j]<0)
1968 gmx_fatal(FARGS,"Found negative temperature in annealing: %f\n",ir->opts.anneal_temp[i][j]);
1972 /* Print out some summary information, to make sure we got it right */
1973 for(i=0,k=0;i<nr;i++) {
1974 if(ir->opts.annealing[i]!=eannNO) {
1975 j = groups->grps[egcTC].nm_ind[i];
1976 fprintf(stderr,"Simulated annealing for group %s: %s, %d timepoints\n",
1977 *(groups->grpname[j]),eann_names[ir->opts.annealing[i]],
1978 ir->opts.anneal_npoints[i]);
1979 fprintf(stderr,"Time (ps) Temperature (K)\n");
1980 /* All terms except the last one */
1981 for(j=0;j<(ir->opts.anneal_npoints[i]-1);j++)
1982 fprintf(stderr,"%9.1f %5.1f\n",ir->opts.anneal_time[i][j],ir->opts.anneal_temp[i][j]);
1984 /* Finally the last one */
1985 j = ir->opts.anneal_npoints[i]-1;
1986 if(ir->opts.annealing[i]==eannSINGLE)
1987 fprintf(stderr,"%9.1f- %5.1f\n",ir->opts.anneal_time[i][j],ir->opts.anneal_temp[i][j]);
1989 fprintf(stderr,"%9.1f %5.1f\n",ir->opts.anneal_time[i][j],ir->opts.anneal_temp[i][j]);
1990 if(fabs(ir->opts.anneal_temp[i][j]-ir->opts.anneal_temp[i][0])>GMX_REAL_EPS)
1991 warning_note(wi,"There is a temperature jump when your annealing loops back.\n");
1999 if (ir->ePull != epullNO) {
2000 make_pull_groups(ir->pull,pull_grp,grps,gnames);
2003 nacc = str_nelem(acc,MAXPTR,ptr1);
2004 nacg = str_nelem(accgrps,MAXPTR,ptr2);
2005 if (nacg*DIM != nacc)
2006 gmx_fatal(FARGS,"Invalid Acceleration input: %d groups and %d acc. values",
2008 do_numbering(natoms,groups,nacg,ptr2,grps,gnames,egcACC,
2009 restnm,egrptpALL_GENREST,bVerbose,wi);
2010 nr = groups->grps[egcACC].nr;
2011 snew(ir->opts.acc,nr);
2014 for(i=k=0; (i<nacg); i++)
2015 for(j=0; (j<DIM); j++,k++)
2016 ir->opts.acc[i][j]=strtod(ptr1[k],NULL);
2018 for(j=0; (j<DIM); j++)
2019 ir->opts.acc[i][j]=0;
2021 nfrdim = str_nelem(frdim,MAXPTR,ptr1);
2022 nfreeze = str_nelem(freeze,MAXPTR,ptr2);
2023 if (nfrdim != DIM*nfreeze)
2024 gmx_fatal(FARGS,"Invalid Freezing input: %d groups and %d freeze values",
2026 do_numbering(natoms,groups,nfreeze,ptr2,grps,gnames,egcFREEZE,
2027 restnm,egrptpALL_GENREST,bVerbose,wi);
2028 nr = groups->grps[egcFREEZE].nr;
2030 snew(ir->opts.nFreeze,nr);
2031 for(i=k=0; (i<nfreeze); i++)
2032 for(j=0; (j<DIM); j++,k++) {
2033 ir->opts.nFreeze[i][j]=(gmx_strncasecmp(ptr1[k],"Y",1)==0);
2034 if (!ir->opts.nFreeze[i][j]) {
2035 if (gmx_strncasecmp(ptr1[k],"N",1) != 0) {
2036 sprintf(warnbuf,"Please use Y(ES) or N(O) for freezedim only "
2037 "(not %s)", ptr1[k]);
2038 warning(wi,warn_buf);
2043 for(j=0; (j<DIM); j++)
2044 ir->opts.nFreeze[i][j]=0;
2046 nenergy=str_nelem(energy,MAXPTR,ptr1);
2047 do_numbering(natoms,groups,nenergy,ptr1,grps,gnames,egcENER,
2048 restnm,egrptpALL_GENREST,bVerbose,wi);
2049 add_wall_energrps(groups,ir->nwall,symtab);
2050 ir->opts.ngener = groups->grps[egcENER].nr;
2051 nvcm=str_nelem(vcm,MAXPTR,ptr1);
2053 do_numbering(natoms,groups,nvcm,ptr1,grps,gnames,egcVCM,
2054 restnm,nvcm==0 ? egrptpALL_GENREST : egrptpPART,bVerbose,wi);
2056 warning(wi,"Some atoms are not part of any center of mass motion removal group.\n"
2057 "This may lead to artifacts.\n"
2058 "In most cases one should use one group for the whole system.");
2061 /* Now we have filled the freeze struct, so we can calculate NRDF */
2062 calc_nrdf(mtop,ir,gnames);
2067 /* Must check per group! */
2068 for(i=0; (i<ir->opts.ngtc); i++)
2069 ntot += ir->opts.nrdf[i];
2070 if (ntot != (DIM*natoms)) {
2071 fac = sqrt(ntot/(DIM*natoms));
2073 fprintf(stderr,"Scaling velocities by a factor of %.3f to account for constraints\n"
2074 "and removal of center of mass motion\n",fac);
2075 for(i=0; (i<natoms); i++)
2076 svmul(fac,v[i],v[i]);
2080 nuser=str_nelem(user1,MAXPTR,ptr1);
2081 do_numbering(natoms,groups,nuser,ptr1,grps,gnames,egcUser1,
2082 restnm,egrptpALL_GENREST,bVerbose,wi);
2083 nuser=str_nelem(user2,MAXPTR,ptr1);
2084 do_numbering(natoms,groups,nuser,ptr1,grps,gnames,egcUser2,
2085 restnm,egrptpALL_GENREST,bVerbose,wi);
2086 nuser=str_nelem(xtc_grps,MAXPTR,ptr1);
2087 do_numbering(natoms,groups,nuser,ptr1,grps,gnames,egcXTC,
2088 restnm,egrptpONE,bVerbose,wi);
2089 nofg = str_nelem(orirefitgrp,MAXPTR,ptr1);
2090 do_numbering(natoms,groups,nofg,ptr1,grps,gnames,egcORFIT,
2091 restnm,egrptpALL_GENREST,bVerbose,wi);
2093 /* QMMM input processing */
2094 nQMg = str_nelem(QMMM,MAXPTR,ptr1);
2095 nQMmethod = str_nelem(QMmethod,MAXPTR,ptr2);
2096 nQMbasis = str_nelem(QMbasis,MAXPTR,ptr3);
2097 if((nQMmethod != nQMg)||(nQMbasis != nQMg)){
2098 gmx_fatal(FARGS,"Invalid QMMM input: %d groups %d basissets"
2099 " and %d methods\n",nQMg,nQMbasis,nQMmethod);
2101 /* group rest, if any, is always MM! */
2102 do_numbering(natoms,groups,nQMg,ptr1,grps,gnames,egcQMMM,
2103 restnm,egrptpALL_GENREST,bVerbose,wi);
2104 nr = nQMg; /*atoms->grps[egcQMMM].nr;*/
2105 ir->opts.ngQM = nQMg;
2106 snew(ir->opts.QMmethod,nr);
2107 snew(ir->opts.QMbasis,nr);
2109 /* input consists of strings: RHF CASSCF PM3 .. These need to be
2110 * converted to the corresponding enum in names.c
2112 ir->opts.QMmethod[i] = search_QMstring(ptr2[i],eQMmethodNR,
2114 ir->opts.QMbasis[i] = search_QMstring(ptr3[i],eQMbasisNR,
2118 nQMmult = str_nelem(QMmult,MAXPTR,ptr1);
2119 nQMcharge = str_nelem(QMcharge,MAXPTR,ptr2);
2120 nbSH = str_nelem(bSH,MAXPTR,ptr3);
2121 snew(ir->opts.QMmult,nr);
2122 snew(ir->opts.QMcharge,nr);
2123 snew(ir->opts.bSH,nr);
2126 ir->opts.QMmult[i] = strtol(ptr1[i],NULL,10);
2127 ir->opts.QMcharge[i] = strtol(ptr2[i],NULL,10);
2128 ir->opts.bSH[i] = (gmx_strncasecmp(ptr3[i],"Y",1)==0);
2131 nCASelec = str_nelem(CASelectrons,MAXPTR,ptr1);
2132 nCASorb = str_nelem(CASorbitals,MAXPTR,ptr2);
2133 snew(ir->opts.CASelectrons,nr);
2134 snew(ir->opts.CASorbitals,nr);
2136 ir->opts.CASelectrons[i]= strtol(ptr1[i],NULL,10);
2137 ir->opts.CASorbitals[i] = strtol(ptr2[i],NULL,10);
2139 /* special optimization options */
2141 nbOPT = str_nelem(bOPT,MAXPTR,ptr1);
2142 nbTS = str_nelem(bTS,MAXPTR,ptr2);
2143 snew(ir->opts.bOPT,nr);
2144 snew(ir->opts.bTS,nr);
2146 ir->opts.bOPT[i] = (gmx_strncasecmp(ptr1[i],"Y",1)==0);
2147 ir->opts.bTS[i] = (gmx_strncasecmp(ptr2[i],"Y",1)==0);
2149 nSAon = str_nelem(SAon,MAXPTR,ptr1);
2150 nSAoff = str_nelem(SAoff,MAXPTR,ptr2);
2151 nSAsteps = str_nelem(SAsteps,MAXPTR,ptr3);
2152 snew(ir->opts.SAon,nr);
2153 snew(ir->opts.SAoff,nr);
2154 snew(ir->opts.SAsteps,nr);
2157 ir->opts.SAon[i] = strtod(ptr1[i],NULL);
2158 ir->opts.SAoff[i] = strtod(ptr2[i],NULL);
2159 ir->opts.SAsteps[i] = strtol(ptr3[i],NULL,10);
2161 /* end of QMMM input */
2164 for(i=0; (i<egcNR); i++) {
2165 fprintf(stderr,"%-16s has %d element(s):",gtypes[i],groups->grps[i].nr);
2166 for(j=0; (j<groups->grps[i].nr); j++)
2167 fprintf(stderr," %s",*(groups->grpname[groups->grps[i].nm_ind[j]]));
2168 fprintf(stderr,"\n");
2171 nr = groups->grps[egcENER].nr;
2172 snew(ir->opts.egp_flags,nr*nr);
2174 bExcl = do_egp_flag(ir,groups,"energygrp-excl",egpexcl,EGP_EXCL);
2175 if (bExcl && EEL_FULL(ir->coulombtype))
2176 warning(wi,"Can not exclude the lattice Coulomb energy between energy groups");
2178 bTable = do_egp_flag(ir,groups,"energygrp-table",egptable,EGP_TABLE);
2179 if (bTable && !(ir->vdwtype == evdwUSER) &&
2180 !(ir->coulombtype == eelUSER) && !(ir->coulombtype == eelPMEUSER) &&
2181 !(ir->coulombtype == eelPMEUSERSWITCH))
2182 gmx_fatal(FARGS,"Can only have energy group pair tables in combination with user tables for VdW and/or Coulomb");
2184 decode_cos(efield_x,&(ir->ex[XX]),FALSE);
2185 decode_cos(efield_xt,&(ir->et[XX]),TRUE);
2186 decode_cos(efield_y,&(ir->ex[YY]),FALSE);
2187 decode_cos(efield_yt,&(ir->et[YY]),TRUE);
2188 decode_cos(efield_z,&(ir->ex[ZZ]),FALSE);
2189 decode_cos(efield_zt,&(ir->et[ZZ]),TRUE);
2191 for(i=0; (i<grps->nr); i++)
2201 static void check_disre(gmx_mtop_t *mtop)
2203 gmx_ffparams_t *ffparams;
2204 t_functype *functype;
2206 int i,ndouble,ftype;
2207 int label,old_label;
2209 if (gmx_mtop_ftype_count(mtop,F_DISRES) > 0) {
2210 ffparams = &mtop->ffparams;
2211 functype = ffparams->functype;
2212 ip = ffparams->iparams;
2215 for(i=0; i<ffparams->ntypes; i++) {
2216 ftype = functype[i];
2217 if (ftype == F_DISRES) {
2218 label = ip[i].disres.label;
2219 if (label == old_label) {
2220 fprintf(stderr,"Distance restraint index %d occurs twice\n",label);
2227 gmx_fatal(FARGS,"Found %d double distance restraint indices,\n"
2228 "probably the parameters for multiple pairs in one restraint "
2229 "are not identical\n",ndouble);
2233 static gmx_bool absolute_reference(t_inputrec *ir,gmx_mtop_t *sys,
2234 gmx_bool posres_only,
2238 gmx_mtop_ilistloop_t iloop;
2248 for(d=0; d<DIM; d++)
2250 AbsRef[d] = (d < ndof_com(ir) ? 0 : 1);
2252 /* Check for freeze groups */
2253 for(g=0; g<ir->opts.ngfrz; g++)
2255 for(d=0; d<DIM; d++)
2257 if (ir->opts.nFreeze[g][d] != 0)
2265 /* Check for position restraints */
2266 iloop = gmx_mtop_ilistloop_init(sys);
2267 while (gmx_mtop_ilistloop_next(iloop,&ilist,&nmol))
2270 (AbsRef[XX] == 0 || AbsRef[YY] == 0 || AbsRef[ZZ] == 0))
2272 for(i=0; i<ilist[F_POSRES].nr; i+=2)
2274 pr = &sys->ffparams.iparams[ilist[F_POSRES].iatoms[i]];
2275 for(d=0; d<DIM; d++)
2277 if (pr->posres.fcA[d] != 0)
2286 return (AbsRef[XX] != 0 && AbsRef[YY] != 0 && AbsRef[ZZ] != 0);
2289 void triple_check(const char *mdparin,t_inputrec *ir,gmx_mtop_t *sys,
2293 int i,m,g,nmol,npct;
2294 gmx_bool bCharge,bAcc;
2295 real gdt_max,*mgrp,mt;
2297 gmx_mtop_atomloop_block_t aloopb;
2298 gmx_mtop_atomloop_all_t aloop;
2301 char warn_buf[STRLEN];
2303 set_warning_line(wi,mdparin,-1);
2305 if (EI_DYNAMICS(ir->eI) && !EI_SD(ir->eI) && ir->eI != eiBD &&
2306 ir->comm_mode == ecmNO &&
2307 !(absolute_reference(ir,sys,FALSE,AbsRef) || ir->nsteps <= 10)) {
2308 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");
2311 /* Check for pressure coupling with absolute position restraints */
2312 if (ir->epc != epcNO && ir->refcoord_scaling == erscNO)
2314 absolute_reference(ir,sys,TRUE,AbsRef);
2316 for(m=0; m<DIM; m++)
2318 if (AbsRef[m] && norm2(ir->compress[m]) > 0)
2320 warning(wi,"You are using pressure coupling with absolute position restraints, this will give artifacts. Use the refcoord_scaling option.");
2328 aloopb = gmx_mtop_atomloop_block_init(sys);
2329 while (gmx_mtop_atomloop_block_next(aloopb,&atom,&nmol)) {
2330 if (atom->q != 0 || atom->qB != 0) {
2336 if (EEL_FULL(ir->coulombtype)) {
2338 "You are using full electrostatics treatment %s for a system without charges.\n"
2339 "This costs a lot of performance for just processing zeros, consider using %s instead.\n",
2340 EELTYPE(ir->coulombtype),EELTYPE(eelCUT));
2341 warning(wi,err_buf);
2344 if (ir->coulombtype == eelCUT && ir->rcoulomb > 0 && !ir->implicit_solvent) {
2346 "You are using a plain Coulomb cut-off, which might produce artifacts.\n"
2347 "You might want to consider using %s electrostatics.\n",
2349 warning_note(wi,err_buf);
2353 /* Generalized reaction field */
2354 if (ir->opts.ngtc == 0) {
2355 sprintf(err_buf,"No temperature coupling while using coulombtype %s",
2357 CHECK(ir->coulombtype == eelGRF);
2360 sprintf(err_buf,"When using coulombtype = %s"
2361 " ref_t for temperature coupling should be > 0",
2363 CHECK((ir->coulombtype == eelGRF) && (ir->opts.ref_t[0] <= 0));
2366 if (ir->eI == eiSD1) {
2368 for(i=0; (i<ir->opts.ngtc); i++)
2369 gdt_max = max(gdt_max,ir->delta_t/ir->opts.tau_t[i]);
2370 if (0.5*gdt_max > 0.0015) {
2371 sprintf(warn_buf,"The relative error with integrator %s is 0.5*delta-t/tau-t = %g, you might want to switch to integrator %s\n",
2372 ei_names[ir->eI],0.5*gdt_max,ei_names[eiSD2]);
2373 warning_note(wi,warn_buf);
2378 for(i=0; (i<sys->groups.grps[egcACC].nr); i++) {
2379 for(m=0; (m<DIM); m++) {
2380 if (fabs(ir->opts.acc[i][m]) > 1e-6) {
2387 snew(mgrp,sys->groups.grps[egcACC].nr);
2388 aloop = gmx_mtop_atomloop_all_init(sys);
2389 while (gmx_mtop_atomloop_all_next(aloop,&i,&atom)) {
2390 mgrp[ggrpnr(&sys->groups,egcACC,i)] += atom->m;
2393 for(i=0; (i<sys->groups.grps[egcACC].nr); i++) {
2394 for(m=0; (m<DIM); m++)
2395 acc[m] += ir->opts.acc[i][m]*mgrp[i];
2398 for(m=0; (m<DIM); m++) {
2399 if (fabs(acc[m]) > 1e-6) {
2400 const char *dim[DIM] = { "X", "Y", "Z" };
2402 "Net Acceleration in %s direction, will %s be corrected\n",
2403 dim[m],ir->nstcomm != 0 ? "" : "not");
2404 if (ir->nstcomm != 0 && m < ndof_com(ir)) {
2406 for (i=0; (i<sys->groups.grps[egcACC].nr); i++)
2407 ir->opts.acc[i][m] -= acc[m];
2414 if (ir->efep != efepNO && ir->sc_alpha != 0 &&
2415 !gmx_within_tol(sys->ffparams.reppow,12.0,10*GMX_DOUBLE_EPS)) {
2416 gmx_fatal(FARGS,"Soft-core interactions are only supported with VdW repulsion power 12");
2419 if (ir->ePull != epullNO) {
2420 if (ir->pull->grp[0].nat == 0) {
2421 absolute_reference(ir,sys,FALSE,AbsRef);
2422 for(m=0; m<DIM; m++) {
2423 if (ir->pull->dim[m] && !AbsRef[m]) {
2424 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.");
2430 if (ir->pull->eGeom == epullgDIRPBC) {
2431 for(i=0; i<3; i++) {
2432 for(m=0; m<=i; m++) {
2433 if ((ir->epc != epcNO && ir->compress[i][m] != 0) ||
2434 ir->deform[i][m] != 0) {
2435 for(g=1; g<ir->pull->ngrp; g++) {
2436 if (ir->pull->grp[g].vec[m] != 0) {
2437 gmx_fatal(FARGS,"Can not have dynamic box while using pull geometry '%s' (dim %c)",EPULLGEOM(ir->pull->eGeom),'x'+m);
2449 void double_check(t_inputrec *ir,matrix box,gmx_bool bConstr,warninp_t wi)
2453 char warn_buf[STRLEN];
2456 ptr = check_box(ir->ePBC,box);
2458 warning_error(wi,ptr);
2461 if (bConstr && ir->eConstrAlg == econtSHAKE) {
2462 if (ir->shake_tol <= 0.0) {
2463 sprintf(warn_buf,"ERROR: shake-tol must be > 0 instead of %g\n",
2465 warning_error(wi,warn_buf);
2468 if (IR_TWINRANGE(*ir) && ir->nstlist > 1) {
2469 sprintf(warn_buf,"With twin-range cut-off's and SHAKE the virial and the pressure are incorrect.");
2470 if (ir->epc == epcNO) {
2471 warning(wi,warn_buf);
2473 warning_error(wi,warn_buf);
2478 if( (ir->eConstrAlg == econtLINCS) && bConstr) {
2479 /* If we have Lincs constraints: */
2480 if(ir->eI==eiMD && ir->etc==etcNO &&
2481 ir->eConstrAlg==econtLINCS && ir->nLincsIter==1) {
2482 sprintf(warn_buf,"For energy conservation with LINCS, lincs_iter should be 2 or larger.\n");
2483 warning_note(wi,warn_buf);
2486 if ((ir->eI == eiCG || ir->eI == eiLBFGS) && (ir->nProjOrder<8)) {
2487 sprintf(warn_buf,"For accurate %s with LINCS constraints, lincs-order should be 8 or more.",ei_names[ir->eI]);
2488 warning_note(wi,warn_buf);
2490 if (ir->epc==epcMTTK) {
2491 warning_error(wi,"MTTK not compatible with lincs -- use shake instead.");
2495 if (ir->LincsWarnAngle > 90.0) {
2496 sprintf(warn_buf,"lincs-warnangle can not be larger than 90 degrees, setting it to 90.\n");
2497 warning(wi,warn_buf);
2498 ir->LincsWarnAngle = 90.0;
2501 if (ir->ePBC != epbcNONE) {
2502 if (ir->nstlist == 0) {
2503 warning(wi,"With nstlist=0 atoms are only put into the box at step 0, therefore drifting atoms might cause the simulation to crash.");
2505 bTWIN = (ir->rlistlong > ir->rlist);
2506 if (ir->ns_type == ensGRID) {
2507 if (sqr(ir->rlistlong) >= max_cutoff2(ir->ePBC,box)) {
2508 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",
2509 bTWIN ? (ir->rcoulomb==ir->rlistlong ? "rcoulomb" : "rvdw"):"rlist");
2510 warning_error(wi,warn_buf);
2513 min_size = min(box[XX][XX],min(box[YY][YY],box[ZZ][ZZ]));
2514 if (2*ir->rlistlong >= min_size) {
2515 sprintf(warn_buf,"ERROR: One of the box lengths is smaller than twice the cut-off length. Increase the box size or decrease rlist.");
2516 warning_error(wi,warn_buf);
2518 fprintf(stderr,"Grid search might allow larger cut-off's than simple search with triclinic boxes.");
2524 void check_chargegroup_radii(const gmx_mtop_t *mtop,const t_inputrec *ir,
2528 real rvdw1,rvdw2,rcoul1,rcoul2;
2529 char warn_buf[STRLEN];
2531 calc_chargegroup_radii(mtop,x,&rvdw1,&rvdw2,&rcoul1,&rcoul2);
2535 printf("Largest charge group radii for Van der Waals: %5.3f, %5.3f nm\n",
2540 printf("Largest charge group radii for Coulomb: %5.3f, %5.3f nm\n",
2546 if (rvdw1 + rvdw2 > ir->rlist ||
2547 rcoul1 + rcoul2 > ir->rlist)
2549 sprintf(warn_buf,"The sum of the two largest charge group radii (%f) is larger than rlist (%f)\n",max(rvdw1+rvdw2,rcoul1+rcoul2),ir->rlist);
2550 warning(wi,warn_buf);
2554 /* Here we do not use the zero at cut-off macro,
2555 * since user defined interactions might purposely
2556 * not be zero at the cut-off.
2558 if (EVDW_IS_ZERO_AT_CUTOFF(ir->vdwtype) &&
2559 rvdw1 + rvdw2 > ir->rlist - ir->rvdw)
2561 sprintf(warn_buf,"The sum of the two largest charge group radii (%f) is larger than rlist (%f) - rvdw (%f)\n",
2563 ir->rlist,ir->rvdw);
2566 warning(wi,warn_buf);
2570 warning_note(wi,warn_buf);
2573 if (EEL_IS_ZERO_AT_CUTOFF(ir->coulombtype) &&
2574 rcoul1 + rcoul2 > ir->rlistlong - ir->rcoulomb)
2576 sprintf(warn_buf,"The sum of the two largest charge group radii (%f) is larger than %s (%f) - rcoulomb (%f)\n",
2578 ir->rlistlong > ir->rlist ? "rlistlong" : "rlist",
2579 ir->rlistlong,ir->rcoulomb);
2582 warning(wi,warn_buf);
2586 warning_note(wi,warn_buf);