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);
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)
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 void check_ir(const char *mdparin,t_inputrec *ir, t_gromppopts *opts,
191 /* Check internal consistency */
193 /* Strange macro: first one fills the err_buf, and then one can check
194 * the condition, which will print the message and increase the error
197 #define CHECK(b) _low_check(b,err_buf,wi)
198 char err_buf[256],warn_buf[STRLEN];
204 t_lambda *fep = ir->fepvals;
205 t_expanded *expand = ir->expandedvals;
207 set_warning_line(wi,mdparin,-1);
209 /* BASIC CUT-OFF STUFF */
210 if (ir->rcoulomb < 0)
212 warning_error(wi,"rcoulomb should be >= 0");
216 warning_error(wi,"rvdw should be >= 0");
220 warning_error(wi,"rlist should be >= 0");
222 if (ir->rlist == 0 ||
223 !((EEL_MIGHT_BE_ZERO_AT_CUTOFF(ir->coulombtype) && ir->rcoulomb > ir->rlist) ||
224 (EVDW_MIGHT_BE_ZERO_AT_CUTOFF(ir->vdwtype) && ir->rvdw > ir->rlist))) {
225 /* No switched potential and/or no twin-range:
226 * we can set the long-range cut-off to the maximum of the other cut-offs.
228 ir->rlistlong = max_cutoff(ir->rlist,max_cutoff(ir->rvdw,ir->rcoulomb));
229 } else if (ir->rlistlong < 0) {
230 ir->rlistlong = max_cutoff(ir->rlist,max_cutoff(ir->rvdw,ir->rcoulomb));
231 sprintf(warn_buf,"rlistlong was not set, setting it to %g (no buffer)",
233 warning(wi,warn_buf);
235 if (ir->rlistlong == 0 && ir->ePBC != epbcNONE) {
236 warning_error(wi,"Can not have an infinite cut-off with PBC");
238 if (ir->rlistlong > 0 && (ir->rlist == 0 || ir->rlistlong < ir->rlist)) {
239 warning_error(wi,"rlistlong can not be shorter than rlist");
241 if (IR_TWINRANGE(*ir) && ir->nstlist <= 0) {
242 warning_error(wi,"Can not have nstlist<=0 with twin-range interactions");
245 /* GENERAL INTEGRATOR STUFF */
246 if (!(ir->eI == eiMD || EI_VV(ir->eI)))
250 if (ir->eI == eiVVAK) {
251 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]);
252 warning_note(wi,warn_buf);
254 if (!EI_DYNAMICS(ir->eI))
258 if (EI_DYNAMICS(ir->eI))
260 if (ir->nstcalcenergy < 0)
262 ir->nstcalcenergy = ir_optimal_nstcalcenergy(ir);
263 if (ir->nstenergy != 0 && ir->nstenergy < ir->nstcalcenergy)
265 /* nstcalcenergy larger than nstener does not make sense.
266 * We ideally want nstcalcenergy=nstener.
270 ir->nstcalcenergy = lcd(ir->nstenergy,ir->nstlist);
274 ir->nstcalcenergy = ir->nstenergy;
278 if (ir->epc != epcNO)
280 if (ir->nstpcouple < 0)
282 ir->nstpcouple = ir_optimal_nstpcouple(ir);
285 if (IR_TWINRANGE(*ir))
287 check_nst("nstlist",ir->nstlist,
288 "nstcalcenergy",&ir->nstcalcenergy,wi);
289 if (ir->epc != epcNO)
291 check_nst("nstlist",ir->nstlist,
292 "nstpcouple",&ir->nstpcouple,wi);
296 if (ir->nstcalcenergy > 1)
298 /* for storing exact averages nstenergy should be
299 * a multiple of nstcalcenergy
301 check_nst("nstcalcenergy",ir->nstcalcenergy,
302 "nstenergy",&ir->nstenergy,wi);
303 if (ir->efep != efepNO)
305 /* nstdhdl should be a multiple of nstcalcenergy */
306 check_nst("nstcalcenergy",ir->nstcalcenergy,
307 "nstdhdl",&ir->fepvals->nstdhdl,wi);
313 if ((EI_SD(ir->eI) || ir->eI == eiBD) &&
314 ir->bContinuation && ir->ld_seed != -1) {
315 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)");
319 if (EI_TPI(ir->eI)) {
320 sprintf(err_buf,"TPI only works with pbc = %s",epbc_names[epbcXYZ]);
321 CHECK(ir->ePBC != epbcXYZ);
322 sprintf(err_buf,"TPI only works with ns = %s",ens_names[ensGRID]);
323 CHECK(ir->ns_type != ensGRID);
324 sprintf(err_buf,"with TPI nstlist should be larger than zero");
325 CHECK(ir->nstlist <= 0);
326 sprintf(err_buf,"TPI does not work with full electrostatics other than PME");
327 CHECK(EEL_FULL(ir->coulombtype) && !EEL_PME(ir->coulombtype));
331 if ( (opts->nshake > 0) && (opts->bMorse) ) {
333 "Using morse bond-potentials while constraining bonds is useless");
334 warning(wi,warn_buf);
337 if ((EI_SD(ir->eI) || ir->eI == eiBD) &&
338 ir->bContinuation && ir->ld_seed != -1) {
339 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)");
341 /* verify simulated tempering options */
344 gmx_bool bAllTempZero = TRUE;
345 for (i=0;i<fep->n_lambda;i++)
347 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]);
348 CHECK((fep->all_lambda[efptTEMPERATURE][i] < 0) || (fep->all_lambda[efptTEMPERATURE][i] > 1));
349 if (fep->all_lambda[efptTEMPERATURE][i] > 0)
351 bAllTempZero = FALSE;
354 sprintf(err_buf,"if simulated tempering is on, temperature-lambdas may not be all zero");
355 CHECK(bAllTempZero==TRUE);
357 sprintf(err_buf,"Simulated tempering is currently only compatible with md-vv");
358 CHECK(ir->eI != eiVV);
360 /* check compatability of the temperature coupling with simulated tempering */
362 if (ir->etc == etcNOSEHOOVER) {
363 sprintf(warn_buf,"Nose-Hoover based temperature control such as [%s] my not be entirelyconsistent with simulated tempering",etcoupl_names[ir->etc]);
364 warning_note(wi,warn_buf);
367 /* check that the temperatures make sense */
369 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);
370 CHECK(ir->simtempvals->simtemp_high <= ir->simtempvals->simtemp_low);
372 sprintf(err_buf,"Higher simulated tempering temperature (%g) must be >= zero",ir->simtempvals->simtemp_high);
373 CHECK(ir->simtempvals->simtemp_high <= 0);
375 sprintf(err_buf,"Lower simulated tempering temperature (%g) must be >= zero",ir->simtempvals->simtemp_low);
376 CHECK(ir->simtempvals->simtemp_low <= 0);
379 /* verify free energy options */
381 if (ir->efep != efepNO) {
383 sprintf(err_buf,"The soft-core power is %d and can only be 1 or 2",
385 CHECK(fep->sc_alpha!=0 && fep->sc_power!=1 && fep->sc_power!=2);
387 sprintf(err_buf,"The soft-core sc-r-power is %d and can only be 6 or 48",
388 (int)fep->sc_r_power);
389 CHECK(fep->sc_alpha!=0 && fep->sc_r_power!=6.0 && fep->sc_r_power!=48.0);
391 /* check validity of options */
392 if (fep->n_lambda > 0 && ir->rlist < max(ir->rvdw,ir->rcoulomb))
395 "For foreign lambda free energy differences it is assumed that the soft-core interactions have no effect beyond the neighborlist cut-off");
396 warning(wi,warn_buf);
399 sprintf(err_buf,"Can't use postive delta-lambda (%g) if initial state/lambda does not start at zero",fep->delta_lambda);
400 CHECK(fep->delta_lambda > 0 && ((fep->init_fep_state !=0) || (fep->init_lambda !=0)));
402 sprintf(err_buf,"Can't use postive delta-lambda (%g) with expanded ensemble simulations",fep->delta_lambda);
403 CHECK(fep->delta_lambda > 0 && (ir->efep == efepEXPANDED));
405 sprintf(err_buf,"Free-energy not implemented for Ewald");
406 CHECK(ir->coulombtype==eelEWALD);
408 /* check validty of lambda inputs */
409 sprintf(err_buf,"initial thermodynamic state %d does not exist, only goes to %d",fep->init_fep_state,fep->n_lambda);
410 CHECK((fep->init_fep_state > fep->n_lambda));
412 for (j=0;j<efptNR;j++)
414 for (i=0;i<fep->n_lambda;i++)
416 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]);
417 CHECK((fep->all_lambda[j][i] < 0) || (fep->all_lambda[j][i] > 1));
421 if ((fep->sc_alpha>0) && (!fep->bScCoul))
423 for (i=0;i<fep->n_lambda;i++)
425 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],
426 fep->all_lambda[efptCOUL][i]);
427 CHECK((fep->sc_alpha>0) &&
428 (((fep->all_lambda[efptCOUL][i] > 0.0) &&
429 (fep->all_lambda[efptCOUL][i] < 1.0)) &&
430 ((fep->all_lambda[efptVDW][i] > 0.0) &&
431 (fep->all_lambda[efptVDW][i] < 1.0))));
435 if ((fep->bScCoul) && (EEL_PME(ir->coulombtype)))
437 sprintf(warn_buf,"With coulomb soft core, the reciprocal space calculation will not necessarily cancel. It may be necessary to decrease the reciprocal space energy, and increase the cutoff radius to get sufficiently close matches to energies with free energy turned off.");
438 warning(wi, warn_buf);
441 /* Free Energy Checks -- In an ideal world, slow growth and FEP would
442 be treated differently, but that's the next step */
444 for (i=0;i<efptNR;i++) {
445 for (j=0;j<fep->n_lambda;j++) {
446 sprintf(err_buf,"%s[%d] must be between 0 and 1",efpt_names[i],j);
447 CHECK((fep->all_lambda[i][j] < 0) || (fep->all_lambda[i][j] > 1));
452 if ((ir->bSimTemp) || (ir->efep == efepEXPANDED)) {
454 expand = ir->expandedvals;
456 /* checking equilibration of weights inputs for validity */
458 sprintf(err_buf,"weight-equil-number-all-lambda (%d) is ignored if lmc-weights-equil is not equal to %s",
459 expand->equil_n_at_lam,elmceq_names[elmceqNUMATLAM]);
460 CHECK((expand->equil_n_at_lam>0) && (expand->elmceq!=elmceqNUMATLAM));
462 sprintf(err_buf,"weight-equil-number-samples (%d) is ignored if lmc-weights-equil is not equal to %s",
463 expand->equil_samples,elmceq_names[elmceqSAMPLES]);
464 CHECK((expand->equil_samples>0) && (expand->elmceq!=elmceqSAMPLES));
466 sprintf(err_buf,"weight-equil-number-steps (%d) is ignored if lmc-weights-equil is not equal to %s",
467 expand->equil_steps,elmceq_names[elmceqSTEPS]);
468 CHECK((expand->equil_steps>0) && (expand->elmceq!=elmceqSTEPS));
470 sprintf(err_buf,"weight-equil-wl-delta (%d) is ignored if lmc-weights-equil is not equal to %s",
471 expand->equil_samples,elmceq_names[elmceqWLDELTA]);
472 CHECK((expand->equil_wl_delta>0) && (expand->elmceq!=elmceqWLDELTA));
474 sprintf(err_buf,"weight-equil-count-ratio (%f) is ignored if lmc-weights-equil is not equal to %s",
475 expand->equil_ratio,elmceq_names[elmceqRATIO]);
476 CHECK((expand->equil_ratio>0) && (expand->elmceq!=elmceqRATIO));
478 sprintf(err_buf,"weight-equil-number-all-lambda (%d) must be a positive integer if lmc-weights-equil=%s",
479 expand->equil_n_at_lam,elmceq_names[elmceqNUMATLAM]);
480 CHECK((expand->equil_n_at_lam<=0) && (expand->elmceq==elmceqNUMATLAM));
482 sprintf(err_buf,"weight-equil-number-samples (%d) must be a positive integer if lmc-weights-equil=%s",
483 expand->equil_samples,elmceq_names[elmceqSAMPLES]);
484 CHECK((expand->equil_samples<=0) && (expand->elmceq==elmceqSAMPLES));
486 sprintf(err_buf,"weight-equil-number-steps (%d) must be a positive integer if lmc-weights-equil=%s",
487 expand->equil_steps,elmceq_names[elmceqSTEPS]);
488 CHECK((expand->equil_steps<=0) && (expand->elmceq==elmceqSTEPS));
490 sprintf(err_buf,"weight-equil-wl-delta (%f) must be > 0 if lmc-weights-equil=%s",
491 expand->equil_wl_delta,elmceq_names[elmceqWLDELTA]);
492 CHECK((expand->equil_wl_delta<=0) && (expand->elmceq==elmceqWLDELTA));
494 sprintf(err_buf,"weight-equil-count-ratio (%f) must be > 0 if lmc-weights-equil=%s",
495 expand->equil_ratio,elmceq_names[elmceqRATIO]);
496 CHECK((expand->equil_ratio<=0) && (expand->elmceq==elmceqRATIO));
498 sprintf(err_buf,"lmc-weights-equil=%s only possible when lmc-stats = %s or lmc-stats %s",
499 elmceq_names[elmceqWLDELTA],elamstats_names[elamstatsWL],elamstats_names[elamstatsWWL]);
500 CHECK((expand->elmceq==elmceqWLDELTA) && (!EWL(expand->elamstats)));
502 sprintf(err_buf,"lmc-repeats (%d) must be greater than 0",expand->lmc_repeats);
503 CHECK((expand->lmc_repeats <= 0));
504 sprintf(err_buf,"minimum-var-min (%d) must be greater than 0",expand->minvarmin);
505 CHECK((expand->minvarmin <= 0));
506 sprintf(err_buf,"weight-c-range (%d) must be greater or equal to 0",expand->c_range);
507 CHECK((expand->c_range < 0));
508 sprintf(err_buf,"init-lambda-state (%d) must be zero if lmc-forced-nstart (%d)> 0 and lmc-move != 'no'",
509 fep->init_fep_state, expand->lmc_forced_nstart);
510 CHECK((fep->init_fep_state!=0) && (expand->lmc_forced_nstart>0) && (expand->elmcmove!=elmcmoveNO));
511 sprintf(err_buf,"lmc-forced-nstart (%d) must not be negative",expand->lmc_forced_nstart);
512 CHECK((expand->lmc_forced_nstart < 0));
513 sprintf(err_buf,"init-lambda-state (%d) must be in the interval [0,number of lambdas)",fep->init_fep_state);
514 CHECK((fep->init_fep_state < 0) || (fep->init_fep_state >= fep->n_lambda));
516 sprintf(err_buf,"init-wl-delta (%f) must be greater than or equal to 0",expand->init_wl_delta);
517 CHECK((expand->init_wl_delta < 0));
518 sprintf(err_buf,"wl-ratio (%f) must be between 0 and 1",expand->wl_ratio);
519 CHECK((expand->wl_ratio <= 0) || (expand->wl_ratio >= 1));
520 sprintf(err_buf,"wl-scale (%f) must be between 0 and 1",expand->wl_scale);
521 CHECK((expand->wl_scale <= 0) || (expand->wl_scale >= 1));
523 /* if there is no temperature control, we need to specify an MC temperature */
524 sprintf(err_buf,"If there is no temperature control, and lmc-mcmove!= 'no',mc_temperature must be set to a positive number");
525 if (expand->nstTij > 0)
527 sprintf(err_buf,"nst-transition-matrix (%d) must be an integer multiple of nstlog (%d)",
528 expand->nstTij,ir->nstlog);
529 CHECK((mod(expand->nstTij,ir->nstlog)!=0));
534 sprintf(err_buf,"walls only work with pbc=%s",epbc_names[epbcXY]);
535 CHECK(ir->nwall && ir->ePBC!=epbcXY);
538 if (ir->ePBC != epbcXYZ && ir->nwall != 2) {
539 if (ir->ePBC == epbcNONE) {
540 if (ir->epc != epcNO) {
541 warning(wi,"Turning off pressure coupling for vacuum system");
545 sprintf(err_buf,"Can not have pressure coupling with pbc=%s",
546 epbc_names[ir->ePBC]);
547 CHECK(ir->epc != epcNO);
549 sprintf(err_buf,"Can not have Ewald with pbc=%s",epbc_names[ir->ePBC]);
550 CHECK(EEL_FULL(ir->coulombtype));
552 sprintf(err_buf,"Can not have dispersion correction with pbc=%s",
553 epbc_names[ir->ePBC]);
554 CHECK(ir->eDispCorr != edispcNO);
557 if (ir->rlist == 0.0) {
558 sprintf(err_buf,"can only have neighborlist cut-off zero (=infinite)\n"
559 "with coulombtype = %s or coulombtype = %s\n"
560 "without periodic boundary conditions (pbc = %s) and\n"
561 "rcoulomb and rvdw set to zero",
562 eel_names[eelCUT],eel_names[eelUSER],epbc_names[epbcNONE]);
563 CHECK(((ir->coulombtype != eelCUT) && (ir->coulombtype != eelUSER)) ||
564 (ir->ePBC != epbcNONE) ||
565 (ir->rcoulomb != 0.0) || (ir->rvdw != 0.0));
567 if (ir->nstlist < 0) {
568 warning_error(wi,"Can not have heuristic neighborlist updates without cut-off");
570 if (ir->nstlist > 0) {
571 warning_note(wi,"Simulating without cut-offs is usually (slightly) faster with nstlist=0, nstype=simple and particle decomposition");
576 if (ir->nstcomm == 0) {
577 ir->comm_mode = ecmNO;
579 if (ir->comm_mode != ecmNO) {
580 if (ir->nstcomm < 0) {
581 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");
582 ir->nstcomm = abs(ir->nstcomm);
585 if (ir->nstcalcenergy > 0 && ir->nstcomm < ir->nstcalcenergy) {
586 warning_note(wi,"nstcomm < nstcalcenergy defeats the purpose of nstcalcenergy, setting nstcomm to nstcalcenergy");
587 ir->nstcomm = ir->nstcalcenergy;
590 if (ir->comm_mode == ecmANGULAR) {
591 sprintf(err_buf,"Can not remove the rotation around the center of mass with periodic molecules");
592 CHECK(ir->bPeriodicMols);
593 if (ir->ePBC != epbcNONE)
594 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).");
598 if (EI_STATE_VELOCITY(ir->eI) && ir->ePBC == epbcNONE && ir->comm_mode != ecmANGULAR) {
599 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.");
602 sprintf(err_buf,"Twin-range neighbour searching (NS) with simple NS"
603 " algorithm not implemented");
604 CHECK(((ir->rcoulomb > ir->rlist) || (ir->rvdw > ir->rlist))
605 && (ir->ns_type == ensSIMPLE));
607 /* TEMPERATURE COUPLING */
608 if (ir->etc == etcYES)
610 ir->etc = etcBERENDSEN;
611 warning_note(wi,"Old option for temperature coupling given: "
612 "changing \"yes\" to \"Berendsen\"\n");
615 if ((ir->etc == etcNOSEHOOVER) || (ir->epc == epcMTTK))
617 if (ir->opts.nhchainlength < 1)
619 sprintf(warn_buf,"number of Nose-Hoover chains (currently %d) cannot be less than 1,reset to 1\n",ir->opts.nhchainlength);
620 ir->opts.nhchainlength =1;
621 warning(wi,warn_buf);
624 if (ir->etc==etcNOSEHOOVER && !EI_VV(ir->eI) && ir->opts.nhchainlength > 1)
626 warning_note(wi,"leapfrog does not yet support Nose-Hoover chains, nhchainlength reset to 1");
627 ir->opts.nhchainlength = 1;
632 ir->opts.nhchainlength = 0;
635 if (ir->eI == eiVVAK) {
636 sprintf(err_buf,"%s implemented primarily for validation, and requires nsttcouple = 1 and nstpcouple = 1.",
638 CHECK((ir->nsttcouple != 1) || (ir->nstpcouple != 1));
641 if (ETC_ANDERSEN(ir->etc))
643 sprintf(err_buf,"%s temperature control not supported for integrator %s.",etcoupl_names[ir->etc],ei_names[ir->eI]);
644 CHECK(!(EI_VV(ir->eI)));
646 for (i=0;i<ir->opts.ngtc;i++)
648 sprintf(err_buf,"all tau_t must currently be equal using Andersen temperature control, violated for group %d",i);
649 CHECK(ir->opts.tau_t[0] != ir->opts.tau_t[i]);
650 sprintf(err_buf,"all tau_t must be postive using Andersen temperature control, tau_t[%d]=%10.6f",
651 i,ir->opts.tau_t[i]);
652 CHECK(ir->opts.tau_t[i]<0);
654 if (ir->nstcomm > 0 && (ir->etc == etcANDERSEN)) {
655 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]);
656 warning_note(wi,warn_buf);
659 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]);
660 CHECK(ir->nstcomm > 1 && (ir->etc == etcANDERSEN));
662 for (i=0;i<ir->opts.ngtc;i++)
664 int nsteps = (int)(ir->opts.tau_t[i]/ir->delta_t);
665 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);
666 CHECK((nsteps % ir->nstcomm) && (ir->etc == etcANDERSENMASSIVE));
669 if (ir->etc == etcBERENDSEN)
671 sprintf(warn_buf,"The %s thermostat does not generate the correct kinetic energy distribution. You might want to consider using the %s thermostat.",
672 ETCOUPLTYPE(ir->etc),ETCOUPLTYPE(etcVRESCALE));
673 warning_note(wi,warn_buf);
676 if ((ir->etc==etcNOSEHOOVER || ETC_ANDERSEN(ir->etc))
677 && ir->epc==epcBERENDSEN)
679 sprintf(warn_buf,"Using Berendsen pressure coupling invalidates the "
680 "true ensemble for the thermostat");
681 warning(wi,warn_buf);
684 /* PRESSURE COUPLING */
685 if (ir->epc == epcISOTROPIC)
687 ir->epc = epcBERENDSEN;
688 warning_note(wi,"Old option for pressure coupling given: "
689 "changing \"Isotropic\" to \"Berendsen\"\n");
692 if (ir->epc != epcNO)
694 dt_pcoupl = ir->nstpcouple*ir->delta_t;
696 sprintf(err_buf,"tau-p must be > 0 instead of %g\n",ir->tau_p);
697 CHECK(ir->tau_p <= 0);
699 if (ir->tau_p/dt_pcoupl < pcouple_min_integration_steps(ir->epc))
701 sprintf(warn_buf,"For proper integration of the %s barostat, tau-p (%g) should be at least %d times larger than nstpcouple*dt (%g)",
702 EPCOUPLTYPE(ir->epc),ir->tau_p,pcouple_min_integration_steps(ir->epc),dt_pcoupl);
703 warning(wi,warn_buf);
706 sprintf(err_buf,"compressibility must be > 0 when using pressure"
707 " coupling %s\n",EPCOUPLTYPE(ir->epc));
708 CHECK(ir->compress[XX][XX] < 0 || ir->compress[YY][YY] < 0 ||
709 ir->compress[ZZ][ZZ] < 0 ||
710 (trace(ir->compress) == 0 && ir->compress[YY][XX] <= 0 &&
711 ir->compress[ZZ][XX] <= 0 && ir->compress[ZZ][YY] <= 0));
713 if (epcPARRINELLORAHMAN == ir->epct && opts->bGenVel)
716 "You are generating velocities so I am assuming you "
717 "are equilibrating a system. You are using "
718 "Parrinello-Rahman pressure coupling, but this can be "
719 "unstable for equilibration. If your system crashes, try "
720 "equilibrating first with Berendsen pressure coupling. If "
721 "you are not equilibrating the system, you can probably "
722 "ignore this warning.");
723 warning(wi,warn_buf);
731 if ((ir->epc!=epcBERENDSEN) && (ir->epc!=epcMTTK))
733 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.");
739 /* More checks are in triple check (grompp.c) */
741 if (ir->coulombtype == eelSWITCH) {
742 sprintf(warn_buf,"coulombtype = %s is only for testing purposes and can lead to serious artifacts, advice: use coulombtype = %s",
743 eel_names[ir->coulombtype],
744 eel_names[eelRF_ZERO]);
745 warning(wi,warn_buf);
748 if (ir->epsilon_r!=1 && ir->implicit_solvent==eisGBSA) {
749 sprintf(warn_buf,"epsilon-r = %g with GB implicit solvent, will use this value for inner dielectric",ir->epsilon_r);
750 warning_note(wi,warn_buf);
753 if (EEL_RF(ir->coulombtype) && ir->epsilon_rf==1 && ir->epsilon_r!=1) {
754 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);
755 warning(wi,warn_buf);
756 ir->epsilon_rf = ir->epsilon_r;
760 if (getenv("GALACTIC_DYNAMICS") == NULL) {
761 sprintf(err_buf,"epsilon-r must be >= 0 instead of %g\n",ir->epsilon_r);
762 CHECK(ir->epsilon_r < 0);
765 if (EEL_RF(ir->coulombtype)) {
766 /* reaction field (at the cut-off) */
768 if (ir->coulombtype == eelRF_ZERO) {
769 sprintf(err_buf,"With coulombtype = %s, epsilon-rf must be 0",
770 eel_names[ir->coulombtype]);
771 CHECK(ir->epsilon_rf != 0);
774 sprintf(err_buf,"epsilon-rf must be >= epsilon-r");
775 CHECK((ir->epsilon_rf < ir->epsilon_r && ir->epsilon_rf != 0) ||
776 (ir->epsilon_r == 0));
777 if (ir->epsilon_rf == ir->epsilon_r) {
778 sprintf(warn_buf,"Using epsilon-rf = epsilon-r with %s does not make sense",
779 eel_names[ir->coulombtype]);
780 warning(wi,warn_buf);
783 /* Allow rlist>rcoulomb for tabulated long range stuff. This just
784 * means the interaction is zero outside rcoulomb, but it helps to
785 * provide accurate energy conservation.
787 if (EEL_MIGHT_BE_ZERO_AT_CUTOFF(ir->coulombtype)) {
788 if (EEL_SWITCHED(ir->coulombtype)) {
790 "With coulombtype = %s rcoulomb_switch must be < rcoulomb",
791 eel_names[ir->coulombtype]);
792 CHECK(ir->rcoulomb_switch >= ir->rcoulomb);
794 } else if (ir->coulombtype == eelCUT || EEL_RF(ir->coulombtype)) {
795 sprintf(err_buf,"With coulombtype = %s, rcoulomb must be >= rlist",
796 eel_names[ir->coulombtype]);
797 CHECK(ir->rlist > ir->rcoulomb);
800 if (EEL_FULL(ir->coulombtype)) {
801 if (ir->coulombtype==eelPMESWITCH || ir->coulombtype==eelPMEUSER ||
802 ir->coulombtype==eelPMEUSERSWITCH) {
803 sprintf(err_buf,"With coulombtype = %s, rcoulomb must be <= rlist",
804 eel_names[ir->coulombtype]);
805 CHECK(ir->rcoulomb > ir->rlist);
807 if (ir->coulombtype == eelPME || ir->coulombtype == eelP3M_AD) {
809 "With coulombtype = %s, rcoulomb must be equal to rlist\n"
810 "If you want optimal energy conservation or exact integration use %s",
811 eel_names[ir->coulombtype],eel_names[eelPMESWITCH]);
814 "With coulombtype = %s, rcoulomb must be equal to rlist",
815 eel_names[ir->coulombtype]);
817 CHECK(ir->rcoulomb != ir->rlist);
821 if (EEL_PME(ir->coulombtype)) {
822 if (ir->pme_order < 3) {
823 warning_error(wi,"pme-order can not be smaller than 3");
827 if (ir->nwall==2 && EEL_FULL(ir->coulombtype)) {
828 if (ir->ewald_geometry == eewg3D) {
829 sprintf(warn_buf,"With pbc=%s you should use ewald-geometry=%s",
830 epbc_names[ir->ePBC],eewg_names[eewg3DC]);
831 warning(wi,warn_buf);
833 /* This check avoids extra pbc coding for exclusion corrections */
834 sprintf(err_buf,"wall-ewald-zfac should be >= 2");
835 CHECK(ir->wall_ewald_zfac < 2);
838 if (EVDW_SWITCHED(ir->vdwtype)) {
839 sprintf(err_buf,"With vdwtype = %s rvdw-switch must be < rvdw",
840 evdw_names[ir->vdwtype]);
841 CHECK(ir->rvdw_switch >= ir->rvdw);
842 } else if (ir->vdwtype == evdwCUT) {
843 sprintf(err_buf,"With vdwtype = %s, rvdw must be >= rlist",evdw_names[ir->vdwtype]);
844 CHECK(ir->rlist > ir->rvdw);
846 if (EEL_IS_ZERO_AT_CUTOFF(ir->coulombtype)
847 && (ir->rlistlong <= ir->rcoulomb)) {
848 sprintf(warn_buf,"For energy conservation with switch/shift potentials, %s should be 0.1 to 0.3 nm larger than rcoulomb.",
849 IR_TWINRANGE(*ir) ? "rlistlong" : "rlist");
850 warning_note(wi,warn_buf);
852 if (EVDW_SWITCHED(ir->vdwtype) && (ir->rlistlong <= ir->rvdw)) {
853 sprintf(warn_buf,"For energy conservation with switch/shift potentials, %s should be 0.1 to 0.3 nm larger than rvdw.",
854 IR_TWINRANGE(*ir) ? "rlistlong" : "rlist");
855 warning_note(wi,warn_buf);
858 if (ir->vdwtype == evdwUSER && ir->eDispCorr != edispcNO) {
859 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.");
862 if (ir->nstlist == -1) {
864 "nstlist=-1 only works with switched or shifted potentials,\n"
865 "suggestion: use vdw-type=%s and coulomb-type=%s",
866 evdw_names[evdwSHIFT],eel_names[eelPMESWITCH]);
867 CHECK(!(EEL_MIGHT_BE_ZERO_AT_CUTOFF(ir->coulombtype) &&
868 EVDW_MIGHT_BE_ZERO_AT_CUTOFF(ir->vdwtype)));
870 sprintf(err_buf,"With nstlist=-1 rvdw and rcoulomb should be smaller than rlist to account for diffusion and possibly charge-group radii");
871 CHECK(ir->rvdw >= ir->rlist || ir->rcoulomb >= ir->rlist);
873 sprintf(err_buf,"nstlist can not be smaller than -1");
874 CHECK(ir->nstlist < -1);
876 if (ir->eI == eiLBFGS && (ir->coulombtype==eelCUT || ir->vdwtype==evdwCUT)
878 warning(wi,"For efficient BFGS minimization, use switch/shift/pme instead of cut-off.");
881 if (ir->eI == eiLBFGS && ir->nbfgscorr <= 0) {
882 warning(wi,"Using L-BFGS with nbfgscorr<=0 just gets you steepest descent.");
885 /* ENERGY CONSERVATION */
888 if (!EVDW_MIGHT_BE_ZERO_AT_CUTOFF(ir->vdwtype) && ir->rvdw > 0)
890 sprintf(warn_buf,"You are using a cut-off for VdW interactions with NVE, for good energy conservation use vdwtype = %s (possibly with DispCorr)",
891 evdw_names[evdwSHIFT]);
892 warning_note(wi,warn_buf);
894 if (!EEL_MIGHT_BE_ZERO_AT_CUTOFF(ir->coulombtype) && ir->rcoulomb > 0)
896 sprintf(warn_buf,"You are using a cut-off for electrostatics with NVE, for good energy conservation use coulombtype = %s or %s",
897 eel_names[eelPMESWITCH],eel_names[eelRF_ZERO]);
898 warning_note(wi,warn_buf);
902 /* IMPLICIT SOLVENT */
903 if(ir->coulombtype==eelGB_NOTUSED)
905 ir->coulombtype=eelCUT;
906 ir->implicit_solvent=eisGBSA;
907 fprintf(stderr,"Note: Old option for generalized born electrostatics given:\n"
908 "Changing coulombtype from \"generalized-born\" to \"cut-off\" and instead\n"
909 "setting implicit-solvent value to \"GBSA\" in input section.\n");
912 if(ir->sa_algorithm==esaSTILL)
914 sprintf(err_buf,"Still SA algorithm not available yet, use %s or %s instead\n",esa_names[esaAPPROX],esa_names[esaNO]);
915 CHECK(ir->sa_algorithm == esaSTILL);
918 if(ir->implicit_solvent==eisGBSA)
920 sprintf(err_buf,"With GBSA implicit solvent, rgbradii must be equal to rlist.");
921 CHECK(ir->rgbradii != ir->rlist);
923 if(ir->coulombtype!=eelCUT)
925 sprintf(err_buf,"With GBSA, coulombtype must be equal to %s\n",eel_names[eelCUT]);
926 CHECK(ir->coulombtype!=eelCUT);
928 if(ir->vdwtype!=evdwCUT)
930 sprintf(err_buf,"With GBSA, vdw-type must be equal to %s\n",evdw_names[evdwCUT]);
931 CHECK(ir->vdwtype!=evdwCUT);
935 sprintf(warn_buf,"Using GBSA with nstgbradii<1, setting nstgbradii=1");
936 warning_note(wi,warn_buf);
939 if(ir->sa_algorithm==esaNO)
941 sprintf(warn_buf,"No SA (non-polar) calculation requested together with GB. Are you sure this is what you want?\n");
942 warning_note(wi,warn_buf);
944 if(ir->sa_surface_tension<0 && ir->sa_algorithm!=esaNO)
946 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");
947 warning_note(wi,warn_buf);
949 if(ir->gb_algorithm==egbSTILL)
951 ir->sa_surface_tension = 0.0049 * CAL2JOULE * 100;
955 ir->sa_surface_tension = 0.0054 * CAL2JOULE * 100;
958 if(ir->sa_surface_tension==0 && ir->sa_algorithm!=esaNO)
960 sprintf(err_buf, "Surface tension set to 0 while SA-calculation requested\n");
961 CHECK(ir->sa_surface_tension==0 && ir->sa_algorithm!=esaNO);
966 if (ir->bAdress && !EI_SD(ir->eI)){
967 warning_error(wi,"AdresS simulation supports only stochastic dynamics");
969 if (ir->bAdress && ir->epc != epcNO){
970 warning_error(wi,"AdresS simulation does not support pressure coupling");
972 if (ir->bAdress && (EEL_PME(ir->coulombtype))){
973 warning_error(wi,"AdresS simulation does not support long-range electrostatics");
978 /* count the number of text elemets separated by whitespace in a string.
979 str = the input string
980 maxptr = the maximum number of allowed elements
981 ptr = the output array of pointers to the first character of each element
982 returns: the number of elements. */
983 int str_nelem(const char *str,int maxptr,char *ptr[])
991 while (*copy != '\0') {
993 gmx_fatal(FARGS,"Too many groups on line: '%s' (max is %d)",
998 while ((*copy != '\0') && !isspace(*copy))
1000 if (*copy != '\0') {
1012 /* interpret a number of doubles from a string and put them in an array,
1013 after allocating space for them.
1014 str = the input string
1015 n = the (pre-allocated) number of doubles read
1016 r = the output array of doubles. */
1017 static void parse_n_real(char *str,int *n,real **r)
1022 *n = str_nelem(str,MAXPTR,ptr);
1025 for(i=0; i<*n; i++) {
1026 (*r)[i] = strtod(ptr[i],NULL);
1030 static void do_fep_params(t_inputrec *ir, char fep_lambda[][STRLEN],char weights[STRLEN]) {
1032 int i,j,max_n_lambda,nweights,nfep[efptNR];
1033 t_lambda *fep = ir->fepvals;
1034 t_expanded *expand = ir->expandedvals;
1035 real **count_fep_lambdas;
1036 gmx_bool bOneLambda = TRUE;
1038 snew(count_fep_lambdas,efptNR);
1040 /* FEP input processing */
1041 /* first, identify the number of lambda values for each type.
1042 All that are nonzero must have the same number */
1044 for (i=0;i<efptNR;i++)
1046 parse_n_real(fep_lambda[i],&(nfep[i]),&(count_fep_lambdas[i]));
1049 /* now, determine the number of components. All must be either zero, or equal. */
1052 for (i=0;i<efptNR;i++)
1054 if (nfep[i] > max_n_lambda) {
1055 max_n_lambda = nfep[i]; /* here's a nonzero one. All of them
1056 must have the same number if its not zero.*/
1061 for (i=0;i<efptNR;i++)
1065 ir->fepvals->separate_dvdl[i] = FALSE;
1067 else if (nfep[i] == max_n_lambda)
1069 if (i!=efptTEMPERATURE) /* we treat this differently -- not really a reason to compute the derivative with
1070 respect to the temperature currently */
1072 ir->fepvals->separate_dvdl[i] = TRUE;
1077 gmx_fatal(FARGS,"Number of lambdas (%d) for FEP type %s not equal to number of other types (%d)",
1078 nfep[i],efpt_names[i],max_n_lambda);
1081 /* we don't print out dhdl if the temperature is changing, since we can't correctly define dhdl in this case */
1082 ir->fepvals->separate_dvdl[efptTEMPERATURE] = FALSE;
1084 /* the number of lambdas is the number we've read in, which is either zero
1085 or the same for all */
1086 fep->n_lambda = max_n_lambda;
1088 /* allocate space for the array of lambda values */
1089 snew(fep->all_lambda,efptNR);
1090 /* if init_lambda is defined, we need to set lambda */
1091 if ((fep->init_lambda > 0) && (fep->n_lambda == 0))
1093 ir->fepvals->separate_dvdl[efptFEP] = TRUE;
1095 /* otherwise allocate the space for all of the lambdas, and transfer the data */
1096 for (i=0;i<efptNR;i++)
1098 snew(fep->all_lambda[i],fep->n_lambda);
1099 if (nfep[i] > 0) /* if it's zero, then the count_fep_lambda arrays
1102 for (j=0;j<fep->n_lambda;j++)
1104 fep->all_lambda[i][j] = (double)count_fep_lambdas[i][j];
1106 sfree(count_fep_lambdas[i]);
1109 sfree(count_fep_lambdas);
1111 /* "fep-vals" is either zero or the full number. If zero, we'll need to define fep-lambdas for internal
1112 bookkeeping -- for now, init_lambda */
1114 if ((nfep[efptFEP] == 0) && (fep->init_lambda >= 0) && (fep->init_lambda <= 1))
1116 for (i=0;i<fep->n_lambda;i++)
1118 fep->all_lambda[efptFEP][i] = fep->init_lambda;
1122 /* check to see if only a single component lambda is defined, and soft core is defined.
1123 In this case, turn on coulomb soft core */
1125 if (max_n_lambda == 0)
1131 for (i=0;i<efptNR;i++)
1133 if ((nfep[i] != 0) && (i!=efptFEP))
1139 if ((bOneLambda) && (fep->sc_alpha > 0))
1141 fep->bScCoul = TRUE;
1144 /* Fill in the others with the efptFEP if they are not explicitly
1145 specified (i.e. nfep[i] == 0). This means if fep is not defined,
1146 they are all zero. */
1148 for (i=0;i<efptNR;i++)
1150 if ((nfep[i] == 0) && (i!=efptFEP))
1152 for (j=0;j<fep->n_lambda;j++)
1154 fep->all_lambda[i][j] = fep->all_lambda[efptFEP][j];
1160 /* make it easier if sc_r_power = 48 by increasing it to the 4th power, to be in the right scale. */
1161 if (fep->sc_r_power == 48)
1163 if (fep->sc_alpha > 0.1)
1165 gmx_fatal(FARGS,"sc_alpha (%f) for sc_r_power = 48 should usually be between 0.001 and 0.004", fep->sc_alpha);
1169 expand = ir->expandedvals;
1170 /* now read in the weights */
1171 parse_n_real(weights,&nweights,&(expand->init_lambda_weights));
1174 expand->bInit_weights = FALSE;
1175 snew(expand->init_lambda_weights,fep->n_lambda); /* initialize to zero */
1177 else if (nweights != fep->n_lambda)
1179 gmx_fatal(FARGS,"Number of weights (%d) is not equal to number of lambda values (%d)",
1180 nweights,fep->n_lambda);
1184 expand->bInit_weights = TRUE;
1186 if ((expand->nstexpanded < 0) && (ir->efep != efepNO)) {
1187 expand->nstexpanded = fep->nstdhdl;
1188 /* if you don't specify nstexpanded when doing expanded ensemble free energy calcs, it is set to nstdhdl */
1190 if ((expand->nstexpanded < 0) && ir->bSimTemp) {
1191 expand->nstexpanded = ir->nstlist;
1192 /* if you don't specify nstexpanded when doing expanded ensemble simulated tempering, it is set to nstlist*/
1197 static void do_simtemp_params(t_inputrec *ir) {
1199 snew(ir->simtempvals->temperatures,ir->fepvals->n_lambda);
1200 GetSimTemps(ir->fepvals->n_lambda,ir->simtempvals,ir->fepvals->all_lambda[efptTEMPERATURE]);
1205 static void do_wall_params(t_inputrec *ir,
1206 char *wall_atomtype, char *wall_density,
1210 char *names[MAXPTR];
1213 opts->wall_atomtype[0] = NULL;
1214 opts->wall_atomtype[1] = NULL;
1216 ir->wall_atomtype[0] = -1;
1217 ir->wall_atomtype[1] = -1;
1218 ir->wall_density[0] = 0;
1219 ir->wall_density[1] = 0;
1223 nstr = str_nelem(wall_atomtype,MAXPTR,names);
1224 if (nstr != ir->nwall)
1226 gmx_fatal(FARGS,"Expected %d elements for wall_atomtype, found %d",
1229 for(i=0; i<ir->nwall; i++)
1231 opts->wall_atomtype[i] = strdup(names[i]);
1234 if (ir->wall_type == ewt93 || ir->wall_type == ewt104) {
1235 nstr = str_nelem(wall_density,MAXPTR,names);
1236 if (nstr != ir->nwall)
1238 gmx_fatal(FARGS,"Expected %d elements for wall-density, found %d",ir->nwall,nstr);
1240 for(i=0; i<ir->nwall; i++)
1242 sscanf(names[i],"%lf",&dbl);
1245 gmx_fatal(FARGS,"wall-density[%d] = %f\n",i,dbl);
1247 ir->wall_density[i] = dbl;
1253 static void add_wall_energrps(gmx_groups_t *groups,int nwall,t_symtab *symtab)
1260 srenew(groups->grpname,groups->ngrpname+nwall);
1261 grps = &(groups->grps[egcENER]);
1262 srenew(grps->nm_ind,grps->nr+nwall);
1263 for(i=0; i<nwall; i++) {
1264 sprintf(str,"wall%d",i);
1265 groups->grpname[groups->ngrpname] = put_symtab(symtab,str);
1266 grps->nm_ind[grps->nr++] = groups->ngrpname++;
1271 void read_expandedparams(int *ninp_p,t_inpfile **inp_p,
1272 t_expanded *expand,warninp_t wi)
1280 /* read expanded ensemble parameters */
1281 CCTYPE ("expanded ensemble variables");
1282 ITYPE ("nstexpanded",expand->nstexpanded,-1);
1283 EETYPE("lmc-stats", expand->elamstats, elamstats_names);
1284 EETYPE("lmc-move", expand->elmcmove, elmcmove_names);
1285 EETYPE("lmc-weights-equil",expand->elmceq,elmceq_names);
1286 ITYPE ("weight-equil-number-all-lambda",expand->equil_n_at_lam,-1);
1287 ITYPE ("weight-equil-number-samples",expand->equil_samples,-1);
1288 ITYPE ("weight-equil-number-steps",expand->equil_steps,-1);
1289 RTYPE ("weight-equil-wl-delta",expand->equil_wl_delta,-1);
1290 RTYPE ("weight-equil-count-ratio",expand->equil_ratio,-1);
1291 CCTYPE("Seed for Monte Carlo in lambda space");
1292 ITYPE ("lmc-seed",expand->lmc_seed,-1);
1293 RTYPE ("mc-temperature",expand->mc_temp,-1);
1294 ITYPE ("lmc-repeats",expand->lmc_repeats,1);
1295 ITYPE ("lmc-gibbsdelta",expand->gibbsdeltalam,-1);
1296 ITYPE ("lmc-forced-nstart",expand->lmc_forced_nstart,0);
1297 EETYPE("symmetrized-transition-matrix", expand->bSymmetrizedTMatrix, yesno_names);
1298 ITYPE("nst-transition-matrix", expand->nstTij, -1);
1299 ITYPE ("mininum-var-min",expand->minvarmin, 100); /*default is reasonable */
1300 ITYPE ("weight-c-range",expand->c_range, 0); /* default is just C=0 */
1301 RTYPE ("wl-scale",expand->wl_scale,0.8);
1302 RTYPE ("wl-ratio",expand->wl_ratio,0.8);
1303 RTYPE ("init-wl-delta",expand->init_wl_delta,1.0);
1304 EETYPE("wl-oneovert",expand->bWLoneovert,yesno_names);
1312 void get_ir(const char *mdparin,const char *mdparout,
1313 t_inputrec *ir,t_gromppopts *opts,
1317 double dumdub[2][6];
1321 char warn_buf[STRLEN];
1322 t_lambda *fep = ir->fepvals;
1323 t_expanded *expand = ir->expandedvals;
1325 inp = read_inpfile(mdparin, &ninp, NULL, wi);
1327 snew(dumstr[0],STRLEN);
1328 snew(dumstr[1],STRLEN);
1330 /* remove the following deprecated commands */
1333 REM_TYPE("domain-decomposition");
1334 REM_TYPE("andersen-seed");
1336 REM_TYPE("dihre-fc");
1337 REM_TYPE("dihre-tau");
1338 REM_TYPE("nstdihreout");
1339 REM_TYPE("nstcheckpoint");
1341 /* replace the following commands with the clearer new versions*/
1342 REPL_TYPE("unconstrained-start","continuation");
1343 REPL_TYPE("foreign-lambda","fep-lambdas");
1345 CCTYPE ("VARIOUS PREPROCESSING OPTIONS");
1346 CTYPE ("Preprocessor information: use cpp syntax.");
1347 CTYPE ("e.g.: -I/home/joe/doe -I/home/mary/roe");
1348 STYPE ("include", opts->include, NULL);
1349 CTYPE ("e.g.: -DPOSRES -DFLEXIBLE (note these variable names are case sensitive)");
1350 STYPE ("define", opts->define, NULL);
1352 CCTYPE ("RUN CONTROL PARAMETERS");
1353 EETYPE("integrator", ir->eI, ei_names);
1354 CTYPE ("Start time and timestep in ps");
1355 RTYPE ("tinit", ir->init_t, 0.0);
1356 RTYPE ("dt", ir->delta_t, 0.001);
1357 STEPTYPE ("nsteps", ir->nsteps, 0);
1358 CTYPE ("For exact run continuation or redoing part of a run");
1359 STEPTYPE ("init-step",ir->init_step, 0);
1360 CTYPE ("Part index is updated automatically on checkpointing (keeps files separate)");
1361 ITYPE ("simulation-part", ir->simulation_part, 1);
1362 CTYPE ("mode for center of mass motion removal");
1363 EETYPE("comm-mode", ir->comm_mode, ecm_names);
1364 CTYPE ("number of steps for center of mass motion removal");
1365 ITYPE ("nstcomm", ir->nstcomm, 10);
1366 CTYPE ("group(s) for center of mass motion removal");
1367 STYPE ("comm-grps", vcm, NULL);
1369 CCTYPE ("LANGEVIN DYNAMICS OPTIONS");
1370 CTYPE ("Friction coefficient (amu/ps) and random seed");
1371 RTYPE ("bd-fric", ir->bd_fric, 0.0);
1372 ITYPE ("ld-seed", ir->ld_seed, 1993);
1375 CCTYPE ("ENERGY MINIMIZATION OPTIONS");
1376 CTYPE ("Force tolerance and initial step-size");
1377 RTYPE ("emtol", ir->em_tol, 10.0);
1378 RTYPE ("emstep", ir->em_stepsize,0.01);
1379 CTYPE ("Max number of iterations in relax-shells");
1380 ITYPE ("niter", ir->niter, 20);
1381 CTYPE ("Step size (ps^2) for minimization of flexible constraints");
1382 RTYPE ("fcstep", ir->fc_stepsize, 0);
1383 CTYPE ("Frequency of steepest descents steps when doing CG");
1384 ITYPE ("nstcgsteep", ir->nstcgsteep, 1000);
1385 ITYPE ("nbfgscorr", ir->nbfgscorr, 10);
1387 CCTYPE ("TEST PARTICLE INSERTION OPTIONS");
1388 RTYPE ("rtpi", ir->rtpi, 0.05);
1390 /* Output options */
1391 CCTYPE ("OUTPUT CONTROL OPTIONS");
1392 CTYPE ("Output frequency for coords (x), velocities (v) and forces (f)");
1393 ITYPE ("nstxout", ir->nstxout, 0);
1394 ITYPE ("nstvout", ir->nstvout, 0);
1395 ITYPE ("nstfout", ir->nstfout, 0);
1396 ir->nstcheckpoint = 1000;
1397 CTYPE ("Output frequency for energies to log file and energy file");
1398 ITYPE ("nstlog", ir->nstlog, 1000);
1399 ITYPE ("nstcalcenergy",ir->nstcalcenergy, -1);
1400 ITYPE ("nstenergy", ir->nstenergy, 100);
1401 CTYPE ("Output frequency and precision for .xtc file");
1402 ITYPE ("nstxtcout", ir->nstxtcout, 0);
1403 RTYPE ("xtc-precision",ir->xtcprec, 1000.0);
1404 CTYPE ("This selects the subset of atoms for the .xtc file. You can");
1405 CTYPE ("select multiple groups. By default all atoms will be written.");
1406 STYPE ("xtc-grps", xtc_grps, NULL);
1407 CTYPE ("Selection of energy groups");
1408 STYPE ("energygrps", energy, NULL);
1410 /* Neighbor searching */
1411 CCTYPE ("NEIGHBORSEARCHING PARAMETERS");
1412 CTYPE ("nblist update frequency");
1413 ITYPE ("nstlist", ir->nstlist, 10);
1414 CTYPE ("ns algorithm (simple or grid)");
1415 EETYPE("ns-type", ir->ns_type, ens_names);
1416 /* set ndelta to the optimal value of 2 */
1418 CTYPE ("Periodic boundary conditions: xyz, no, xy");
1419 EETYPE("pbc", ir->ePBC, epbc_names);
1420 EETYPE("periodic-molecules", ir->bPeriodicMols, yesno_names);
1421 CTYPE ("nblist cut-off");
1422 RTYPE ("rlist", ir->rlist, -1);
1423 CTYPE ("long-range cut-off for switched potentials");
1424 RTYPE ("rlistlong", ir->rlistlong, -1);
1426 /* Electrostatics */
1427 CCTYPE ("OPTIONS FOR ELECTROSTATICS AND VDW");
1428 CTYPE ("Method for doing electrostatics");
1429 EETYPE("coulombtype", ir->coulombtype, eel_names);
1430 CTYPE ("cut-off lengths");
1431 RTYPE ("rcoulomb-switch", ir->rcoulomb_switch, 0.0);
1432 RTYPE ("rcoulomb", ir->rcoulomb, -1);
1433 CTYPE ("Relative dielectric constant for the medium and the reaction field");
1434 RTYPE ("epsilon-r", ir->epsilon_r, 1.0);
1435 RTYPE ("epsilon-rf", ir->epsilon_rf, 0.0);
1436 CTYPE ("Method for doing Van der Waals");
1437 EETYPE("vdw-type", ir->vdwtype, evdw_names);
1438 CTYPE ("cut-off lengths");
1439 RTYPE ("rvdw-switch", ir->rvdw_switch, 0.0);
1440 RTYPE ("rvdw", ir->rvdw, -1);
1441 CTYPE ("Apply long range dispersion corrections for Energy and Pressure");
1442 EETYPE("DispCorr", ir->eDispCorr, edispc_names);
1443 CTYPE ("Extension of the potential lookup tables beyond the cut-off");
1444 RTYPE ("table-extension", ir->tabext, 1.0);
1445 CTYPE ("Seperate tables between energy group pairs");
1446 STYPE ("energygrp-table", egptable, NULL);
1447 CTYPE ("Spacing for the PME/PPPM FFT grid");
1448 RTYPE ("fourierspacing", opts->fourierspacing,0.12);
1449 CTYPE ("FFT grid size, when a value is 0 fourierspacing will be used");
1450 ITYPE ("fourier-nx", ir->nkx, 0);
1451 ITYPE ("fourier-ny", ir->nky, 0);
1452 ITYPE ("fourier-nz", ir->nkz, 0);
1453 CTYPE ("EWALD/PME/PPPM parameters");
1454 ITYPE ("pme-order", ir->pme_order, 4);
1455 RTYPE ("ewald-rtol", ir->ewald_rtol, 0.00001);
1456 EETYPE("ewald-geometry", ir->ewald_geometry, eewg_names);
1457 RTYPE ("epsilon-surface", ir->epsilon_surface, 0.0);
1458 EETYPE("optimize-fft",ir->bOptFFT, yesno_names);
1460 CCTYPE("IMPLICIT SOLVENT ALGORITHM");
1461 EETYPE("implicit-solvent", ir->implicit_solvent, eis_names);
1463 CCTYPE ("GENERALIZED BORN ELECTROSTATICS");
1464 CTYPE ("Algorithm for calculating Born radii");
1465 EETYPE("gb-algorithm", ir->gb_algorithm, egb_names);
1466 CTYPE ("Frequency of calculating the Born radii inside rlist");
1467 ITYPE ("nstgbradii", ir->nstgbradii, 1);
1468 CTYPE ("Cutoff for Born radii calculation; the contribution from atoms");
1469 CTYPE ("between rlist and rgbradii is updated every nstlist steps");
1470 RTYPE ("rgbradii", ir->rgbradii, 1.0);
1471 CTYPE ("Dielectric coefficient of the implicit solvent");
1472 RTYPE ("gb-epsilon-solvent",ir->gb_epsilon_solvent, 80.0);
1473 CTYPE ("Salt concentration in M for Generalized Born models");
1474 RTYPE ("gb-saltconc", ir->gb_saltconc, 0.0);
1475 CTYPE ("Scaling factors used in the OBC GB model. Default values are OBC(II)");
1476 RTYPE ("gb-obc-alpha", ir->gb_obc_alpha, 1.0);
1477 RTYPE ("gb-obc-beta", ir->gb_obc_beta, 0.8);
1478 RTYPE ("gb-obc-gamma", ir->gb_obc_gamma, 4.85);
1479 RTYPE ("gb-dielectric-offset", ir->gb_dielectric_offset, 0.009);
1480 EETYPE("sa-algorithm", ir->sa_algorithm, esa_names);
1481 CTYPE ("Surface tension (kJ/mol/nm^2) for the SA (nonpolar surface) part of GBSA");
1482 CTYPE ("The value -1 will set default value for Still/HCT/OBC GB-models.");
1483 RTYPE ("sa-surface-tension", ir->sa_surface_tension, -1);
1485 /* Coupling stuff */
1486 CCTYPE ("OPTIONS FOR WEAK COUPLING ALGORITHMS");
1487 CTYPE ("Temperature coupling");
1488 EETYPE("tcoupl", ir->etc, etcoupl_names);
1489 ITYPE ("nsttcouple", ir->nsttcouple, -1);
1490 ITYPE("nh-chain-length", ir->opts.nhchainlength, NHCHAINLENGTH);
1491 EETYPE("print-nose-hoover-chain-variables", ir->bPrintNHChains, yesno_names);
1492 CTYPE ("Groups to couple separately");
1493 STYPE ("tc-grps", tcgrps, NULL);
1494 CTYPE ("Time constant (ps) and reference temperature (K)");
1495 STYPE ("tau-t", tau_t, NULL);
1496 STYPE ("ref-t", ref_t, NULL);
1497 CTYPE ("pressure coupling");
1498 EETYPE("pcoupl", ir->epc, epcoupl_names);
1499 EETYPE("pcoupltype", ir->epct, epcoupltype_names);
1500 ITYPE ("nstpcouple", ir->nstpcouple, -1);
1501 CTYPE ("Time constant (ps), compressibility (1/bar) and reference P (bar)");
1502 RTYPE ("tau-p", ir->tau_p, 1.0);
1503 STYPE ("compressibility", dumstr[0], NULL);
1504 STYPE ("ref-p", dumstr[1], NULL);
1505 CTYPE ("Scaling of reference coordinates, No, All or COM");
1506 EETYPE ("refcoord-scaling",ir->refcoord_scaling,erefscaling_names);
1509 CCTYPE ("OPTIONS FOR QMMM calculations");
1510 EETYPE("QMMM", ir->bQMMM, yesno_names);
1511 CTYPE ("Groups treated Quantum Mechanically");
1512 STYPE ("QMMM-grps", QMMM, NULL);
1513 CTYPE ("QM method");
1514 STYPE("QMmethod", QMmethod, NULL);
1515 CTYPE ("QMMM scheme");
1516 EETYPE("QMMMscheme", ir->QMMMscheme, eQMMMscheme_names);
1517 CTYPE ("QM basisset");
1518 STYPE("QMbasis", QMbasis, NULL);
1519 CTYPE ("QM charge");
1520 STYPE ("QMcharge", QMcharge,NULL);
1521 CTYPE ("QM multiplicity");
1522 STYPE ("QMmult", QMmult,NULL);
1523 CTYPE ("Surface Hopping");
1524 STYPE ("SH", bSH, NULL);
1525 CTYPE ("CAS space options");
1526 STYPE ("CASorbitals", CASorbitals, NULL);
1527 STYPE ("CASelectrons", CASelectrons, NULL);
1528 STYPE ("SAon", SAon, NULL);
1529 STYPE ("SAoff",SAoff,NULL);
1530 STYPE ("SAsteps", SAsteps, NULL);
1531 CTYPE ("Scale factor for MM charges");
1532 RTYPE ("MMChargeScaleFactor", ir->scalefactor, 1.0);
1533 CTYPE ("Optimization of QM subsystem");
1534 STYPE ("bOPT", bOPT, NULL);
1535 STYPE ("bTS", bTS, NULL);
1537 /* Simulated annealing */
1538 CCTYPE("SIMULATED ANNEALING");
1539 CTYPE ("Type of annealing for each temperature group (no/single/periodic)");
1540 STYPE ("annealing", anneal, NULL);
1541 CTYPE ("Number of time points to use for specifying annealing in each group");
1542 STYPE ("annealing-npoints", anneal_npoints, NULL);
1543 CTYPE ("List of times at the annealing points for each group");
1544 STYPE ("annealing-time", anneal_time, NULL);
1545 CTYPE ("Temp. at each annealing point, for each group.");
1546 STYPE ("annealing-temp", anneal_temp, NULL);
1549 CCTYPE ("GENERATE VELOCITIES FOR STARTUP RUN");
1550 EETYPE("gen-vel", opts->bGenVel, yesno_names);
1551 RTYPE ("gen-temp", opts->tempi, 300.0);
1552 ITYPE ("gen-seed", opts->seed, 173529);
1555 CCTYPE ("OPTIONS FOR BONDS");
1556 EETYPE("constraints", opts->nshake, constraints);
1557 CTYPE ("Type of constraint algorithm");
1558 EETYPE("constraint-algorithm", ir->eConstrAlg, econstr_names);
1559 CTYPE ("Do not constrain the start configuration");
1560 EETYPE("continuation", ir->bContinuation, yesno_names);
1561 CTYPE ("Use successive overrelaxation to reduce the number of shake iterations");
1562 EETYPE("Shake-SOR", ir->bShakeSOR, yesno_names);
1563 CTYPE ("Relative tolerance of shake");
1564 RTYPE ("shake-tol", ir->shake_tol, 0.0001);
1565 CTYPE ("Highest order in the expansion of the constraint coupling matrix");
1566 ITYPE ("lincs-order", ir->nProjOrder, 4);
1567 CTYPE ("Number of iterations in the final step of LINCS. 1 is fine for");
1568 CTYPE ("normal simulations, but use 2 to conserve energy in NVE runs.");
1569 CTYPE ("For energy minimization with constraints it should be 4 to 8.");
1570 ITYPE ("lincs-iter", ir->nLincsIter, 1);
1571 CTYPE ("Lincs will write a warning to the stderr if in one step a bond");
1572 CTYPE ("rotates over more degrees than");
1573 RTYPE ("lincs-warnangle", ir->LincsWarnAngle, 30.0);
1574 CTYPE ("Convert harmonic bonds to morse potentials");
1575 EETYPE("morse", opts->bMorse,yesno_names);
1577 /* Energy group exclusions */
1578 CCTYPE ("ENERGY GROUP EXCLUSIONS");
1579 CTYPE ("Pairs of energy groups for which all non-bonded interactions are excluded");
1580 STYPE ("energygrp-excl", egpexcl, NULL);
1584 CTYPE ("Number of walls, type, atom types, densities and box-z scale factor for Ewald");
1585 ITYPE ("nwall", ir->nwall, 0);
1586 EETYPE("wall-type", ir->wall_type, ewt_names);
1587 RTYPE ("wall-r-linpot", ir->wall_r_linpot, -1);
1588 STYPE ("wall-atomtype", wall_atomtype, NULL);
1589 STYPE ("wall-density", wall_density, NULL);
1590 RTYPE ("wall-ewald-zfac", ir->wall_ewald_zfac, 3);
1593 CCTYPE("COM PULLING");
1594 CTYPE("Pull type: no, umbrella, constraint or constant-force");
1595 EETYPE("pull", ir->ePull, epull_names);
1596 if (ir->ePull != epullNO) {
1598 pull_grp = read_pullparams(&ninp,&inp,ir->pull,&opts->pull_start,wi);
1601 /* Enforced rotation */
1602 CCTYPE("ENFORCED ROTATION");
1603 CTYPE("Enforced rotation: No or Yes");
1604 EETYPE("rotation", ir->bRot, yesno_names);
1607 rot_grp = read_rotparams(&ninp,&inp,ir->rot,wi);
1611 CCTYPE("NMR refinement stuff");
1612 CTYPE ("Distance restraints type: No, Simple or Ensemble");
1613 EETYPE("disre", ir->eDisre, edisre_names);
1614 CTYPE ("Force weighting of pairs in one distance restraint: Conservative or Equal");
1615 EETYPE("disre-weighting", ir->eDisreWeighting, edisreweighting_names);
1616 CTYPE ("Use sqrt of the time averaged times the instantaneous violation");
1617 EETYPE("disre-mixed", ir->bDisreMixed, yesno_names);
1618 RTYPE ("disre-fc", ir->dr_fc, 1000.0);
1619 RTYPE ("disre-tau", ir->dr_tau, 0.0);
1620 CTYPE ("Output frequency for pair distances to energy file");
1621 ITYPE ("nstdisreout", ir->nstdisreout, 100);
1622 CTYPE ("Orientation restraints: No or Yes");
1623 EETYPE("orire", opts->bOrire, yesno_names);
1624 CTYPE ("Orientation restraints force constant and tau for time averaging");
1625 RTYPE ("orire-fc", ir->orires_fc, 0.0);
1626 RTYPE ("orire-tau", ir->orires_tau, 0.0);
1627 STYPE ("orire-fitgrp",orirefitgrp, NULL);
1628 CTYPE ("Output frequency for trace(SD) and S to energy file");
1629 ITYPE ("nstorireout", ir->nstorireout, 100);
1631 /* free energy variables */
1632 CCTYPE ("Free energy variables");
1633 EETYPE("free-energy", ir->efep, efep_names);
1634 STYPE ("couple-moltype", couple_moltype, NULL);
1635 EETYPE("couple-lambda0", opts->couple_lam0, couple_lam);
1636 EETYPE("couple-lambda1", opts->couple_lam1, couple_lam);
1637 EETYPE("couple-intramol", opts->bCoupleIntra, yesno_names);
1639 RTYPE ("init-lambda", fep->init_lambda,-1); /* start with -1 so
1641 it was not entered */
1642 ITYPE ("init-lambda-state", fep->init_fep_state,0);
1643 RTYPE ("delta-lambda",fep->delta_lambda,0.0);
1644 ITYPE ("nstdhdl",fep->nstdhdl, 10);
1645 STYPE ("fep-lambdas", fep_lambda[efptFEP], NULL);
1646 STYPE ("mass-lambdas", fep_lambda[efptMASS], NULL);
1647 STYPE ("coul-lambdas", fep_lambda[efptCOUL], NULL);
1648 STYPE ("vdw-lambdas", fep_lambda[efptVDW], NULL);
1649 STYPE ("bonded-lambdas", fep_lambda[efptBONDED], NULL);
1650 STYPE ("restraint-lambdas", fep_lambda[efptRESTRAINT], NULL);
1651 STYPE ("temperature-lambdas", fep_lambda[efptTEMPERATURE], NULL);
1652 STYPE ("init-lambda-weights",lambda_weights,NULL);
1653 EETYPE("dhdl-print-energy", fep->bPrintEnergy, yesno_names);
1654 RTYPE ("sc-alpha",fep->sc_alpha,0.0);
1655 ITYPE ("sc-power",fep->sc_power,1);
1656 RTYPE ("sc-r-power",fep->sc_r_power,6.0);
1657 RTYPE ("sc-sigma",fep->sc_sigma,0.3);
1658 EETYPE("sc-coul",fep->bScCoul,yesno_names);
1659 ITYPE ("dh_hist_size", fep->dh_hist_size, 0);
1660 RTYPE ("dh_hist_spacing", fep->dh_hist_spacing, 0.1);
1661 EETYPE("separate-dhdl-file", fep->separate_dhdl_file,
1662 separate_dhdl_file_names);
1663 EETYPE("dhdl-derivatives", fep->dhdl_derivatives, dhdl_derivatives_names);
1664 ITYPE ("dh_hist_size", fep->dh_hist_size, 0);
1665 RTYPE ("dh_hist_spacing", fep->dh_hist_spacing, 0.1);
1667 /* Non-equilibrium MD stuff */
1668 CCTYPE("Non-equilibrium MD stuff");
1669 STYPE ("acc-grps", accgrps, NULL);
1670 STYPE ("accelerate", acc, NULL);
1671 STYPE ("freezegrps", freeze, NULL);
1672 STYPE ("freezedim", frdim, NULL);
1673 RTYPE ("cos-acceleration", ir->cos_accel, 0);
1674 STYPE ("deform", deform, NULL);
1676 /* simulated tempering variables */
1677 CCTYPE("simulated tempering variables");
1678 EETYPE("simulated-tempering",ir->bSimTemp,yesno_names);
1679 EETYPE("simulated-tempering-scaling",ir->simtempvals->eSimTempScale,esimtemp_names);
1680 RTYPE("sim-temp-low",ir->simtempvals->simtemp_low,300.0);
1681 RTYPE("sim-temp-high",ir->simtempvals->simtemp_high,300.0);
1683 /* expanded ensemble variables */
1684 if (ir->efep==efepEXPANDED || ir->bSimTemp)
1686 read_expandedparams(&ninp,&inp,expand,wi);
1689 /* Electric fields */
1690 CCTYPE("Electric fields");
1691 CTYPE ("Format is number of terms (int) and for all terms an amplitude (real)");
1692 CTYPE ("and a phase angle (real)");
1693 STYPE ("E-x", efield_x, NULL);
1694 STYPE ("E-xt", efield_xt, NULL);
1695 STYPE ("E-y", efield_y, NULL);
1696 STYPE ("E-yt", efield_yt, NULL);
1697 STYPE ("E-z", efield_z, NULL);
1698 STYPE ("E-zt", efield_zt, NULL);
1700 /* AdResS defined thingies */
1701 CCTYPE ("AdResS parameters");
1702 EETYPE("adress", ir->bAdress, yesno_names);
1705 read_adressparams(&ninp,&inp,ir->adress,wi);
1708 /* User defined thingies */
1709 CCTYPE ("User defined thingies");
1710 STYPE ("user1-grps", user1, NULL);
1711 STYPE ("user2-grps", user2, NULL);
1712 ITYPE ("userint1", ir->userint1, 0);
1713 ITYPE ("userint2", ir->userint2, 0);
1714 ITYPE ("userint3", ir->userint3, 0);
1715 ITYPE ("userint4", ir->userint4, 0);
1716 RTYPE ("userreal1", ir->userreal1, 0);
1717 RTYPE ("userreal2", ir->userreal2, 0);
1718 RTYPE ("userreal3", ir->userreal3, 0);
1719 RTYPE ("userreal4", ir->userreal4, 0);
1722 write_inpfile(mdparout,ninp,inp,FALSE,wi);
1723 for (i=0; (i<ninp); i++) {
1725 sfree(inp[i].value);
1729 /* Process options if necessary */
1730 for(m=0; m<2; m++) {
1731 for(i=0; i<2*DIM; i++)
1736 if (sscanf(dumstr[m],"%lf",&(dumdub[m][XX]))!=1) {
1737 warning_error(wi,"Pressure coupling not enough values (I need 1)");
1739 dumdub[m][YY]=dumdub[m][ZZ]=dumdub[m][XX];
1741 case epctSEMIISOTROPIC:
1742 case epctSURFACETENSION:
1743 if (sscanf(dumstr[m],"%lf%lf",
1744 &(dumdub[m][XX]),&(dumdub[m][ZZ]))!=2) {
1745 warning_error(wi,"Pressure coupling not enough values (I need 2)");
1747 dumdub[m][YY]=dumdub[m][XX];
1749 case epctANISOTROPIC:
1750 if (sscanf(dumstr[m],"%lf%lf%lf%lf%lf%lf",
1751 &(dumdub[m][XX]),&(dumdub[m][YY]),&(dumdub[m][ZZ]),
1752 &(dumdub[m][3]),&(dumdub[m][4]),&(dumdub[m][5]))!=6) {
1753 warning_error(wi,"Pressure coupling not enough values (I need 6)");
1757 gmx_fatal(FARGS,"Pressure coupling type %s not implemented yet",
1758 epcoupltype_names[ir->epct]);
1762 clear_mat(ir->ref_p);
1763 clear_mat(ir->compress);
1764 for(i=0; i<DIM; i++) {
1765 ir->ref_p[i][i] = dumdub[1][i];
1766 ir->compress[i][i] = dumdub[0][i];
1768 if (ir->epct == epctANISOTROPIC) {
1769 ir->ref_p[XX][YY] = dumdub[1][3];
1770 ir->ref_p[XX][ZZ] = dumdub[1][4];
1771 ir->ref_p[YY][ZZ] = dumdub[1][5];
1772 if (ir->ref_p[XX][YY]!=0 && ir->ref_p[XX][ZZ]!=0 && ir->ref_p[YY][ZZ]!=0) {
1773 warning(wi,"All off-diagonal reference pressures are non-zero. Are you sure you want to apply a threefold shear stress?\n");
1775 ir->compress[XX][YY] = dumdub[0][3];
1776 ir->compress[XX][ZZ] = dumdub[0][4];
1777 ir->compress[YY][ZZ] = dumdub[0][5];
1778 for(i=0; i<DIM; i++) {
1779 for(m=0; m<i; m++) {
1780 ir->ref_p[i][m] = ir->ref_p[m][i];
1781 ir->compress[i][m] = ir->compress[m][i];
1786 if (ir->comm_mode == ecmNO)
1789 opts->couple_moltype = NULL;
1790 if (strlen(couple_moltype) > 0)
1792 if (ir->efep != efepNO)
1794 opts->couple_moltype = strdup(couple_moltype);
1795 if (opts->couple_lam0 == opts->couple_lam1)
1797 warning(wi,"The lambda=0 and lambda=1 states for coupling are identical");
1799 if (ir->eI == eiMD && (opts->couple_lam0 == ecouplamNONE ||
1800 opts->couple_lam1 == ecouplamNONE))
1802 warning(wi,"For proper sampling of the (nearly) decoupled state, stochastic dynamics should be used");
1807 warning(wi,"Can not couple a molecule with free_energy = no");
1810 /* FREE ENERGY AND EXPANDED ENSEMBLE OPTIONS */
1811 if (ir->efep != efepNO) {
1812 if (fep->delta_lambda > 0) {
1813 ir->efep = efepSLOWGROWTH;
1818 fep->bPrintEnergy = TRUE;
1819 /* always print out the energy to dhdl if we are doing expanded ensemble, since we need the total energy
1820 if the temperature is changing. */
1823 if ((ir->efep != efepNO) || ir->bSimTemp)
1825 ir->bExpanded = FALSE;
1826 if ((ir->efep == efepEXPANDED) || ir->bSimTemp)
1828 ir->bExpanded = TRUE;
1830 do_fep_params(ir,fep_lambda,lambda_weights);
1831 if (ir->bSimTemp) { /* done after fep params */
1832 do_simtemp_params(ir);
1837 ir->fepvals->n_lambda = 0;
1840 /* WALL PARAMETERS */
1842 do_wall_params(ir,wall_atomtype,wall_density,opts);
1844 /* ORIENTATION RESTRAINT PARAMETERS */
1846 if (opts->bOrire && str_nelem(orirefitgrp,MAXPTR,NULL)!=1) {
1847 warning_error(wi,"ERROR: Need one orientation restraint fit group\n");
1850 /* DEFORMATION PARAMETERS */
1852 clear_mat(ir->deform);
1857 m = sscanf(deform,"%lf %lf %lf %lf %lf %lf",
1858 &(dumdub[0][0]),&(dumdub[0][1]),&(dumdub[0][2]),
1859 &(dumdub[0][3]),&(dumdub[0][4]),&(dumdub[0][5]));
1862 ir->deform[i][i] = dumdub[0][i];
1864 ir->deform[YY][XX] = dumdub[0][3];
1865 ir->deform[ZZ][XX] = dumdub[0][4];
1866 ir->deform[ZZ][YY] = dumdub[0][5];
1867 if (ir->epc != epcNO) {
1870 if (ir->deform[i][j]!=0 && ir->compress[i][j]!=0) {
1871 warning_error(wi,"A box element has deform set and compressibility > 0");
1875 if (ir->deform[i][j]!=0) {
1876 for(m=j; m<DIM; m++)
1877 if (ir->compress[m][j]!=0) {
1878 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.");
1879 warning(wi,warn_buf);
1888 static int search_QMstring(char *s,int ng,const char *gn[])
1890 /* same as normal search_string, but this one searches QM strings */
1893 for(i=0; (i<ng); i++)
1894 if (gmx_strcasecmp(s,gn[i]) == 0)
1897 gmx_fatal(FARGS,"this QM method or basisset (%s) is not implemented\n!",s);
1901 } /* search_QMstring */
1904 int search_string(char *s,int ng,char *gn[])
1908 for(i=0; (i<ng); i++)
1910 if (gmx_strcasecmp(s,gn[i]) == 0)
1917 "Group %s referenced in the .mdp file was not found in the index file.\n"
1918 "Group names must match either [moleculetype] names or custom index group\n"
1919 "names, in which case you must supply an index file to the '-n' option\n"
1926 static gmx_bool do_numbering(int natoms,gmx_groups_t *groups,int ng,char *ptrs[],
1927 t_blocka *block,char *gnames[],
1928 int gtype,int restnm,
1929 int grptp,gmx_bool bVerbose,
1932 unsigned short *cbuf;
1933 t_grps *grps=&(groups->grps[gtype]);
1934 int i,j,gid,aj,ognr,ntot=0;
1937 char warn_buf[STRLEN];
1941 fprintf(debug,"Starting numbering %d groups of type %d\n",ng,gtype);
1944 title = gtypes[gtype];
1947 /* Mark all id's as not set */
1948 for(i=0; (i<natoms); i++)
1953 snew(grps->nm_ind,ng+1); /* +1 for possible rest group */
1954 for(i=0; (i<ng); i++)
1956 /* Lookup the group name in the block structure */
1957 gid = search_string(ptrs[i],block->nr,gnames);
1958 if ((grptp != egrptpONE) || (i == 0))
1960 grps->nm_ind[grps->nr++]=gid;
1964 fprintf(debug,"Found gid %d for group %s\n",gid,ptrs[i]);
1967 /* Now go over the atoms in the group */
1968 for(j=block->index[gid]; (j<block->index[gid+1]); j++)
1973 /* Range checking */
1974 if ((aj < 0) || (aj >= natoms))
1976 gmx_fatal(FARGS,"Invalid atom number %d in indexfile",aj);
1978 /* Lookup up the old group number */
1982 gmx_fatal(FARGS,"Atom %d in multiple %s groups (%d and %d)",
1983 aj+1,title,ognr+1,i+1);
1987 /* Store the group number in buffer */
1988 if (grptp == egrptpONE)
2001 /* Now check whether we have done all atoms */
2005 if (grptp == egrptpALL)
2007 gmx_fatal(FARGS,"%d atoms are not part of any of the %s groups",
2010 else if (grptp == egrptpPART)
2012 sprintf(warn_buf,"%d atoms are not part of any of the %s groups",
2014 warning_note(wi,warn_buf);
2016 /* Assign all atoms currently unassigned to a rest group */
2017 for(j=0; (j<natoms); j++)
2019 if (cbuf[j] == NOGID)
2025 if (grptp != egrptpPART)
2030 "Making dummy/rest group for %s containing %d elements\n",
2033 /* Add group name "rest" */
2034 grps->nm_ind[grps->nr] = restnm;
2036 /* Assign the rest name to all atoms not currently assigned to a group */
2037 for(j=0; (j<natoms); j++)
2039 if (cbuf[j] == NOGID)
2048 if (grps->nr == 1 && (ntot == 0 || ntot == natoms))
2050 /* All atoms are part of one (or no) group, no index required */
2051 groups->ngrpnr[gtype] = 0;
2052 groups->grpnr[gtype] = NULL;
2056 groups->ngrpnr[gtype] = natoms;
2057 snew(groups->grpnr[gtype],natoms);
2058 for(j=0; (j<natoms); j++)
2060 groups->grpnr[gtype][j] = cbuf[j];
2066 return (bRest && grptp == egrptpPART);
2069 static void calc_nrdf(gmx_mtop_t *mtop,t_inputrec *ir,char **gnames)
2072 gmx_groups_t *groups;
2074 int natoms,ai,aj,i,j,d,g,imin,jmin,nc;
2076 int *nrdf2,*na_vcm,na_tot;
2077 double *nrdf_tc,*nrdf_vcm,nrdf_uc,n_sub=0;
2078 gmx_mtop_atomloop_all_t aloop;
2080 int mb,mol,ftype,as;
2081 gmx_molblock_t *molb;
2082 gmx_moltype_t *molt;
2085 * First calc 3xnr-atoms for each group
2086 * then subtract half a degree of freedom for each constraint
2088 * Only atoms and nuclei contribute to the degrees of freedom...
2093 groups = &mtop->groups;
2094 natoms = mtop->natoms;
2096 /* Allocate one more for a possible rest group */
2097 /* We need to sum degrees of freedom into doubles,
2098 * since floats give too low nrdf's above 3 million atoms.
2100 snew(nrdf_tc,groups->grps[egcTC].nr+1);
2101 snew(nrdf_vcm,groups->grps[egcVCM].nr+1);
2102 snew(na_vcm,groups->grps[egcVCM].nr+1);
2104 for(i=0; i<groups->grps[egcTC].nr; i++)
2106 for(i=0; i<groups->grps[egcVCM].nr+1; i++)
2110 aloop = gmx_mtop_atomloop_all_init(mtop);
2111 while (gmx_mtop_atomloop_all_next(aloop,&i,&atom)) {
2113 if (atom->ptype == eptAtom || atom->ptype == eptNucleus) {
2114 g = ggrpnr(groups,egcFREEZE,i);
2115 /* Double count nrdf for particle i */
2116 for(d=0; d<DIM; d++) {
2117 if (opts->nFreeze[g][d] == 0) {
2121 nrdf_tc [ggrpnr(groups,egcTC ,i)] += 0.5*nrdf2[i];
2122 nrdf_vcm[ggrpnr(groups,egcVCM,i)] += 0.5*nrdf2[i];
2127 for(mb=0; mb<mtop->nmolblock; mb++) {
2128 molb = &mtop->molblock[mb];
2129 molt = &mtop->moltype[molb->type];
2130 atom = molt->atoms.atom;
2131 for(mol=0; mol<molb->nmol; mol++) {
2132 for (ftype=F_CONSTR; ftype<=F_CONSTRNC; ftype++) {
2133 ia = molt->ilist[ftype].iatoms;
2134 for(i=0; i<molt->ilist[ftype].nr; ) {
2135 /* Subtract degrees of freedom for the constraints,
2136 * if the particles still have degrees of freedom left.
2137 * If one of the particles is a vsite or a shell, then all
2138 * constraint motion will go there, but since they do not
2139 * contribute to the constraints the degrees of freedom do not
2144 if (((atom[ia[1]].ptype == eptNucleus) ||
2145 (atom[ia[1]].ptype == eptAtom)) &&
2146 ((atom[ia[2]].ptype == eptNucleus) ||
2147 (atom[ia[2]].ptype == eptAtom))) {
2156 imin = min(imin,nrdf2[ai]);
2157 jmin = min(jmin,nrdf2[aj]);
2160 nrdf_tc [ggrpnr(groups,egcTC ,ai)] -= 0.5*imin;
2161 nrdf_tc [ggrpnr(groups,egcTC ,aj)] -= 0.5*jmin;
2162 nrdf_vcm[ggrpnr(groups,egcVCM,ai)] -= 0.5*imin;
2163 nrdf_vcm[ggrpnr(groups,egcVCM,aj)] -= 0.5*jmin;
2165 ia += interaction_function[ftype].nratoms+1;
2166 i += interaction_function[ftype].nratoms+1;
2169 ia = molt->ilist[F_SETTLE].iatoms;
2170 for(i=0; i<molt->ilist[F_SETTLE].nr; ) {
2171 /* Subtract 1 dof from every atom in the SETTLE */
2172 for(j=0; j<3; j++) {
2174 imin = min(2,nrdf2[ai]);
2176 nrdf_tc [ggrpnr(groups,egcTC ,ai)] -= 0.5*imin;
2177 nrdf_vcm[ggrpnr(groups,egcVCM,ai)] -= 0.5*imin;
2182 as += molt->atoms.nr;
2186 if (ir->ePull == epullCONSTRAINT) {
2187 /* Correct nrdf for the COM constraints.
2188 * We correct using the TC and VCM group of the first atom
2189 * in the reference and pull group. If atoms in one pull group
2190 * belong to different TC or VCM groups it is anyhow difficult
2191 * to determine the optimal nrdf assignment.
2194 if (pull->eGeom == epullgPOS) {
2196 for(i=0; i<DIM; i++) {
2203 for(i=0; i<pull->ngrp; i++) {
2205 if (pull->grp[0].nat > 0) {
2206 /* Subtract 1/2 dof from the reference group */
2207 ai = pull->grp[0].ind[0];
2208 if (nrdf_tc[ggrpnr(groups,egcTC,ai)] > 1) {
2209 nrdf_tc [ggrpnr(groups,egcTC ,ai)] -= 0.5;
2210 nrdf_vcm[ggrpnr(groups,egcVCM,ai)] -= 0.5;
2214 /* Subtract 1/2 dof from the pulled group */
2215 ai = pull->grp[1+i].ind[0];
2216 nrdf_tc [ggrpnr(groups,egcTC ,ai)] -= 0.5*imin;
2217 nrdf_vcm[ggrpnr(groups,egcVCM,ai)] -= 0.5*imin;
2218 if (nrdf_tc[ggrpnr(groups,egcTC,ai)] < 0)
2219 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)]]);
2223 if (ir->nstcomm != 0) {
2224 /* Subtract 3 from the number of degrees of freedom in each vcm group
2225 * when com translation is removed and 6 when rotation is removed
2228 switch (ir->comm_mode) {
2230 n_sub = ndof_com(ir);
2237 gmx_incons("Checking comm_mode");
2240 for(i=0; i<groups->grps[egcTC].nr; i++) {
2241 /* Count the number of atoms of TC group i for every VCM group */
2242 for(j=0; j<groups->grps[egcVCM].nr+1; j++)
2245 for(ai=0; ai<natoms; ai++)
2246 if (ggrpnr(groups,egcTC,ai) == i) {
2247 na_vcm[ggrpnr(groups,egcVCM,ai)]++;
2250 /* Correct for VCM removal according to the fraction of each VCM
2251 * group present in this TC group.
2253 nrdf_uc = nrdf_tc[i];
2255 fprintf(debug,"T-group[%d] nrdf_uc = %g, n_sub = %g\n",
2259 for(j=0; j<groups->grps[egcVCM].nr+1; j++) {
2260 if (nrdf_vcm[j] > n_sub) {
2261 nrdf_tc[i] += nrdf_uc*((double)na_vcm[j]/(double)na_tot)*
2262 (nrdf_vcm[j] - n_sub)/nrdf_vcm[j];
2265 fprintf(debug," nrdf_vcm[%d] = %g, nrdf = %g\n",
2266 j,nrdf_vcm[j],nrdf_tc[i]);
2271 for(i=0; (i<groups->grps[egcTC].nr); i++) {
2272 opts->nrdf[i] = nrdf_tc[i];
2273 if (opts->nrdf[i] < 0)
2276 "Number of degrees of freedom in T-Coupling group %s is %.2f\n",
2277 gnames[groups->grps[egcTC].nm_ind[i]],opts->nrdf[i]);
2286 static void decode_cos(char *s,t_cosines *cosine,gmx_bool bTime)
2289 char format[STRLEN],f1[STRLEN];
2300 sscanf(t,"%d",&(cosine->n));
2301 if (cosine->n <= 0) {
2304 snew(cosine->a,cosine->n);
2305 snew(cosine->phi,cosine->n);
2307 sprintf(format,"%%*d");
2308 for(i=0; (i<cosine->n); i++) {
2310 strcat(f1,"%lf%lf");
2311 if (sscanf(t,f1,&a,&phi) < 2)
2312 gmx_fatal(FARGS,"Invalid input for electric field shift: '%s'",t);
2315 strcat(format,"%*lf%*lf");
2322 static gmx_bool do_egp_flag(t_inputrec *ir,gmx_groups_t *groups,
2323 const char *option,const char *val,int flag)
2325 /* The maximum number of energy group pairs would be MAXPTR*(MAXPTR+1)/2.
2326 * But since this is much larger than STRLEN, such a line can not be parsed.
2327 * The real maximum is the number of names that fit in a string: STRLEN/2.
2329 #define EGP_MAX (STRLEN/2)
2331 char *names[EGP_MAX];
2335 gnames = groups->grpname;
2337 nelem = str_nelem(val,EGP_MAX,names);
2339 gmx_fatal(FARGS,"The number of groups for %s is odd",option);
2340 nr = groups->grps[egcENER].nr;
2342 for(i=0; i<nelem/2; i++) {
2345 gmx_strcasecmp(names[2*i],*(gnames[groups->grps[egcENER].nm_ind[j]])))
2348 gmx_fatal(FARGS,"%s in %s is not an energy group\n",
2352 gmx_strcasecmp(names[2*i+1],*(gnames[groups->grps[egcENER].nm_ind[k]])))
2355 gmx_fatal(FARGS,"%s in %s is not an energy group\n",
2356 names[2*i+1],option);
2357 if ((j < nr) && (k < nr)) {
2358 ir->opts.egp_flags[nr*j+k] |= flag;
2359 ir->opts.egp_flags[nr*k+j] |= flag;
2367 void do_index(const char* mdparin, const char *ndx,
2370 t_inputrec *ir,rvec *v,
2374 gmx_groups_t *groups;
2378 char warnbuf[STRLEN],**gnames;
2379 int nr,ntcg,ntau_t,nref_t,nacc,nofg,nSA,nSA_points,nSA_time,nSA_temp;
2382 int nacg,nfreeze,nfrdim,nenergy,nvcm,nuser;
2383 char *ptr1[MAXPTR],*ptr2[MAXPTR],*ptr3[MAXPTR];
2386 gmx_bool bExcl,bTable,bSetTCpar,bAnneal,bRest;
2387 int nQMmethod,nQMbasis,nQMcharge,nQMmult,nbSH,nCASorb,nCASelec,
2388 nSAon,nSAoff,nSAsteps,nQMg,nbOPT,nbTS;
2389 char warn_buf[STRLEN];
2392 fprintf(stderr,"processing index file...\n");
2396 snew(grps->index,1);
2398 atoms_all = gmx_mtop_global_atoms(mtop);
2399 analyse(&atoms_all,grps,&gnames,FALSE,TRUE);
2400 free_t_atoms(&atoms_all,FALSE);
2402 grps = init_index(ndx,&gnames);
2405 groups = &mtop->groups;
2406 natoms = mtop->natoms;
2407 symtab = &mtop->symtab;
2409 snew(groups->grpname,grps->nr+1);
2411 for(i=0; (i<grps->nr); i++) {
2412 groups->grpname[i] = put_symtab(symtab,gnames[i]);
2414 groups->grpname[i] = put_symtab(symtab,"rest");
2416 srenew(gnames,grps->nr+1);
2417 gnames[restnm] = *(groups->grpname[i]);
2418 groups->ngrpname = grps->nr+1;
2420 set_warning_line(wi,mdparin,-1);
2422 ntau_t = str_nelem(tau_t,MAXPTR,ptr1);
2423 nref_t = str_nelem(ref_t,MAXPTR,ptr2);
2424 ntcg = str_nelem(tcgrps,MAXPTR,ptr3);
2425 if ((ntau_t != ntcg) || (nref_t != ntcg)) {
2426 gmx_fatal(FARGS,"Invalid T coupling input: %d groups, %d ref-t values and "
2427 "%d tau-t values",ntcg,nref_t,ntau_t);
2430 bSetTCpar = (ir->etc || EI_SD(ir->eI) || ir->eI==eiBD || EI_TPI(ir->eI));
2431 do_numbering(natoms,groups,ntcg,ptr3,grps,gnames,egcTC,
2432 restnm,bSetTCpar ? egrptpALL : egrptpALL_GENREST,bVerbose,wi);
2433 nr = groups->grps[egcTC].nr;
2435 snew(ir->opts.nrdf,nr);
2436 snew(ir->opts.tau_t,nr);
2437 snew(ir->opts.ref_t,nr);
2438 if (ir->eI==eiBD && ir->bd_fric==0) {
2439 fprintf(stderr,"bd-fric=0, so tau-t will be used as the inverse friction constant(s)\n");
2446 gmx_fatal(FARGS,"Not enough ref-t and tau-t values!");
2450 for(i=0; (i<nr); i++)
2452 ir->opts.tau_t[i] = strtod(ptr1[i],NULL);
2453 if ((ir->eI == eiBD || ir->eI == eiSD2) && ir->opts.tau_t[i] <= 0)
2455 sprintf(warn_buf,"With integrator %s tau-t should be larger than 0",ei_names[ir->eI]);
2456 warning_error(wi,warn_buf);
2458 if ((ir->etc == etcVRESCALE && ir->opts.tau_t[i] >= 0) ||
2459 (ir->etc != etcVRESCALE && ir->opts.tau_t[i] > 0))
2461 tau_min = min(tau_min,ir->opts.tau_t[i]);
2464 if (ir->etc != etcNO && ir->nsttcouple == -1)
2466 ir->nsttcouple = ir_optimal_nsttcouple(ir);
2471 if ((ir->etc==etcNOSEHOOVER) && (ir->epc==epcBERENDSEN)) {
2472 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");
2474 if ((ir->epc==epcMTTK) && (ir->etc>etcNO))
2477 mincouple = ir->nsttcouple;
2478 if (ir->nstpcouple < mincouple)
2480 mincouple = ir->nstpcouple;
2482 ir->nstpcouple = mincouple;
2483 ir->nsttcouple = mincouple;
2484 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);
2485 warning_note(wi,warn_buf);
2488 /* velocity verlet with averaged kinetic energy KE = 0.5*(v(t+1/2) - v(t-1/2)) is implemented
2489 primarily for testing purposes, and does not work with temperature coupling other than 1 */
2491 if (ETC_ANDERSEN(ir->etc)) {
2492 if (ir->nsttcouple != 1) {
2494 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");
2495 warning_note(wi,warn_buf);
2498 nstcmin = tcouple_min_integration_steps(ir->etc);
2501 if (tau_min/(ir->delta_t*ir->nsttcouple) < nstcmin)
2503 sprintf(warn_buf,"For proper integration of the %s thermostat, tau-t (%g) should be at least %d times larger than nsttcouple*dt (%g)",
2504 ETCOUPLTYPE(ir->etc),
2506 ir->nsttcouple*ir->delta_t);
2507 warning(wi,warn_buf);
2510 for(i=0; (i<nr); i++)
2512 ir->opts.ref_t[i] = strtod(ptr2[i],NULL);
2513 if (ir->opts.ref_t[i] < 0)
2515 gmx_fatal(FARGS,"ref-t for group %d negative",i);
2518 /* set the lambda mc temperature to the md integrator temperature (which should be defined
2519 if we are in this conditional) if mc_temp is negative */
2520 if (ir->expandedvals->mc_temp < 0)
2522 ir->expandedvals->mc_temp = ir->opts.ref_t[0]; /*for now, set to the first reft */
2526 /* Simulated annealing for each group. There are nr groups */
2527 nSA = str_nelem(anneal,MAXPTR,ptr1);
2528 if (nSA == 1 && (ptr1[0][0]=='n' || ptr1[0][0]=='N'))
2530 if(nSA>0 && nSA != nr)
2531 gmx_fatal(FARGS,"Not enough annealing values: %d (for %d groups)\n",nSA,nr);
2533 snew(ir->opts.annealing,nr);
2534 snew(ir->opts.anneal_npoints,nr);
2535 snew(ir->opts.anneal_time,nr);
2536 snew(ir->opts.anneal_temp,nr);
2538 ir->opts.annealing[i]=eannNO;
2539 ir->opts.anneal_npoints[i]=0;
2540 ir->opts.anneal_time[i]=NULL;
2541 ir->opts.anneal_temp[i]=NULL;
2546 if(ptr1[i][0]=='n' || ptr1[i][0]=='N') {
2547 ir->opts.annealing[i]=eannNO;
2548 } else if(ptr1[i][0]=='s'|| ptr1[i][0]=='S') {
2549 ir->opts.annealing[i]=eannSINGLE;
2551 } else if(ptr1[i][0]=='p'|| ptr1[i][0]=='P') {
2552 ir->opts.annealing[i]=eannPERIODIC;
2557 /* Read the other fields too */
2558 nSA_points = str_nelem(anneal_npoints,MAXPTR,ptr1);
2560 gmx_fatal(FARGS,"Found %d annealing-npoints values for %d groups\n",nSA_points,nSA);
2561 for(k=0,i=0;i<nr;i++) {
2562 ir->opts.anneal_npoints[i]=strtol(ptr1[i],NULL,10);
2563 if(ir->opts.anneal_npoints[i]==1)
2564 gmx_fatal(FARGS,"Please specify at least a start and an end point for annealing\n");
2565 snew(ir->opts.anneal_time[i],ir->opts.anneal_npoints[i]);
2566 snew(ir->opts.anneal_temp[i],ir->opts.anneal_npoints[i]);
2567 k += ir->opts.anneal_npoints[i];
2570 nSA_time = str_nelem(anneal_time,MAXPTR,ptr1);
2572 gmx_fatal(FARGS,"Found %d annealing-time values, wanter %d\n",nSA_time,k);
2573 nSA_temp = str_nelem(anneal_temp,MAXPTR,ptr2);
2575 gmx_fatal(FARGS,"Found %d annealing-temp values, wanted %d\n",nSA_temp,k);
2577 for(i=0,k=0;i<nr;i++) {
2579 for(j=0;j<ir->opts.anneal_npoints[i];j++) {
2580 ir->opts.anneal_time[i][j]=strtod(ptr1[k],NULL);
2581 ir->opts.anneal_temp[i][j]=strtod(ptr2[k],NULL);
2583 if(ir->opts.anneal_time[i][0] > (ir->init_t+GMX_REAL_EPS))
2584 gmx_fatal(FARGS,"First time point for annealing > init_t.\n");
2587 if(ir->opts.anneal_time[i][j]<ir->opts.anneal_time[i][j-1])
2588 gmx_fatal(FARGS,"Annealing timepoints out of order: t=%f comes after t=%f\n",
2589 ir->opts.anneal_time[i][j],ir->opts.anneal_time[i][j-1]);
2591 if(ir->opts.anneal_temp[i][j]<0)
2592 gmx_fatal(FARGS,"Found negative temperature in annealing: %f\n",ir->opts.anneal_temp[i][j]);
2596 /* Print out some summary information, to make sure we got it right */
2597 for(i=0,k=0;i<nr;i++) {
2598 if(ir->opts.annealing[i]!=eannNO) {
2599 j = groups->grps[egcTC].nm_ind[i];
2600 fprintf(stderr,"Simulated annealing for group %s: %s, %d timepoints\n",
2601 *(groups->grpname[j]),eann_names[ir->opts.annealing[i]],
2602 ir->opts.anneal_npoints[i]);
2603 fprintf(stderr,"Time (ps) Temperature (K)\n");
2604 /* All terms except the last one */
2605 for(j=0;j<(ir->opts.anneal_npoints[i]-1);j++)
2606 fprintf(stderr,"%9.1f %5.1f\n",ir->opts.anneal_time[i][j],ir->opts.anneal_temp[i][j]);
2608 /* Finally the last one */
2609 j = ir->opts.anneal_npoints[i]-1;
2610 if(ir->opts.annealing[i]==eannSINGLE)
2611 fprintf(stderr,"%9.1f- %5.1f\n",ir->opts.anneal_time[i][j],ir->opts.anneal_temp[i][j]);
2613 fprintf(stderr,"%9.1f %5.1f\n",ir->opts.anneal_time[i][j],ir->opts.anneal_temp[i][j]);
2614 if(fabs(ir->opts.anneal_temp[i][j]-ir->opts.anneal_temp[i][0])>GMX_REAL_EPS)
2615 warning_note(wi,"There is a temperature jump when your annealing loops back.\n");
2623 if (ir->ePull != epullNO) {
2624 make_pull_groups(ir->pull,pull_grp,grps,gnames);
2628 make_rotation_groups(ir->rot,rot_grp,grps,gnames);
2631 nacc = str_nelem(acc,MAXPTR,ptr1);
2632 nacg = str_nelem(accgrps,MAXPTR,ptr2);
2633 if (nacg*DIM != nacc)
2634 gmx_fatal(FARGS,"Invalid Acceleration input: %d groups and %d acc. values",
2636 do_numbering(natoms,groups,nacg,ptr2,grps,gnames,egcACC,
2637 restnm,egrptpALL_GENREST,bVerbose,wi);
2638 nr = groups->grps[egcACC].nr;
2639 snew(ir->opts.acc,nr);
2642 for(i=k=0; (i<nacg); i++)
2643 for(j=0; (j<DIM); j++,k++)
2644 ir->opts.acc[i][j]=strtod(ptr1[k],NULL);
2646 for(j=0; (j<DIM); j++)
2647 ir->opts.acc[i][j]=0;
2649 nfrdim = str_nelem(frdim,MAXPTR,ptr1);
2650 nfreeze = str_nelem(freeze,MAXPTR,ptr2);
2651 if (nfrdim != DIM*nfreeze)
2652 gmx_fatal(FARGS,"Invalid Freezing input: %d groups and %d freeze values",
2654 do_numbering(natoms,groups,nfreeze,ptr2,grps,gnames,egcFREEZE,
2655 restnm,egrptpALL_GENREST,bVerbose,wi);
2656 nr = groups->grps[egcFREEZE].nr;
2658 snew(ir->opts.nFreeze,nr);
2659 for(i=k=0; (i<nfreeze); i++)
2660 for(j=0; (j<DIM); j++,k++) {
2661 ir->opts.nFreeze[i][j]=(gmx_strncasecmp(ptr1[k],"Y",1)==0);
2662 if (!ir->opts.nFreeze[i][j]) {
2663 if (gmx_strncasecmp(ptr1[k],"N",1) != 0) {
2664 sprintf(warnbuf,"Please use Y(ES) or N(O) for freezedim only "
2665 "(not %s)", ptr1[k]);
2666 warning(wi,warn_buf);
2671 for(j=0; (j<DIM); j++)
2672 ir->opts.nFreeze[i][j]=0;
2674 nenergy=str_nelem(energy,MAXPTR,ptr1);
2675 do_numbering(natoms,groups,nenergy,ptr1,grps,gnames,egcENER,
2676 restnm,egrptpALL_GENREST,bVerbose,wi);
2677 add_wall_energrps(groups,ir->nwall,symtab);
2678 ir->opts.ngener = groups->grps[egcENER].nr;
2679 nvcm=str_nelem(vcm,MAXPTR,ptr1);
2681 do_numbering(natoms,groups,nvcm,ptr1,grps,gnames,egcVCM,
2682 restnm,nvcm==0 ? egrptpALL_GENREST : egrptpPART,bVerbose,wi);
2684 warning(wi,"Some atoms are not part of any center of mass motion removal group.\n"
2685 "This may lead to artifacts.\n"
2686 "In most cases one should use one group for the whole system.");
2689 /* Now we have filled the freeze struct, so we can calculate NRDF */
2690 calc_nrdf(mtop,ir,gnames);
2695 /* Must check per group! */
2696 for(i=0; (i<ir->opts.ngtc); i++)
2697 ntot += ir->opts.nrdf[i];
2698 if (ntot != (DIM*natoms)) {
2699 fac = sqrt(ntot/(DIM*natoms));
2701 fprintf(stderr,"Scaling velocities by a factor of %.3f to account for constraints\n"
2702 "and removal of center of mass motion\n",fac);
2703 for(i=0; (i<natoms); i++)
2704 svmul(fac,v[i],v[i]);
2708 nuser=str_nelem(user1,MAXPTR,ptr1);
2709 do_numbering(natoms,groups,nuser,ptr1,grps,gnames,egcUser1,
2710 restnm,egrptpALL_GENREST,bVerbose,wi);
2711 nuser=str_nelem(user2,MAXPTR,ptr1);
2712 do_numbering(natoms,groups,nuser,ptr1,grps,gnames,egcUser2,
2713 restnm,egrptpALL_GENREST,bVerbose,wi);
2714 nuser=str_nelem(xtc_grps,MAXPTR,ptr1);
2715 do_numbering(natoms,groups,nuser,ptr1,grps,gnames,egcXTC,
2716 restnm,egrptpONE,bVerbose,wi);
2717 nofg = str_nelem(orirefitgrp,MAXPTR,ptr1);
2718 do_numbering(natoms,groups,nofg,ptr1,grps,gnames,egcORFIT,
2719 restnm,egrptpALL_GENREST,bVerbose,wi);
2721 /* QMMM input processing */
2722 nQMg = str_nelem(QMMM,MAXPTR,ptr1);
2723 nQMmethod = str_nelem(QMmethod,MAXPTR,ptr2);
2724 nQMbasis = str_nelem(QMbasis,MAXPTR,ptr3);
2725 if((nQMmethod != nQMg)||(nQMbasis != nQMg)){
2726 gmx_fatal(FARGS,"Invalid QMMM input: %d groups %d basissets"
2727 " and %d methods\n",nQMg,nQMbasis,nQMmethod);
2729 /* group rest, if any, is always MM! */
2730 do_numbering(natoms,groups,nQMg,ptr1,grps,gnames,egcQMMM,
2731 restnm,egrptpALL_GENREST,bVerbose,wi);
2732 nr = nQMg; /*atoms->grps[egcQMMM].nr;*/
2733 ir->opts.ngQM = nQMg;
2734 snew(ir->opts.QMmethod,nr);
2735 snew(ir->opts.QMbasis,nr);
2737 /* input consists of strings: RHF CASSCF PM3 .. These need to be
2738 * converted to the corresponding enum in names.c
2740 ir->opts.QMmethod[i] = search_QMstring(ptr2[i],eQMmethodNR,
2742 ir->opts.QMbasis[i] = search_QMstring(ptr3[i],eQMbasisNR,
2746 nQMmult = str_nelem(QMmult,MAXPTR,ptr1);
2747 nQMcharge = str_nelem(QMcharge,MAXPTR,ptr2);
2748 nbSH = str_nelem(bSH,MAXPTR,ptr3);
2749 snew(ir->opts.QMmult,nr);
2750 snew(ir->opts.QMcharge,nr);
2751 snew(ir->opts.bSH,nr);
2754 ir->opts.QMmult[i] = strtol(ptr1[i],NULL,10);
2755 ir->opts.QMcharge[i] = strtol(ptr2[i],NULL,10);
2756 ir->opts.bSH[i] = (gmx_strncasecmp(ptr3[i],"Y",1)==0);
2759 nCASelec = str_nelem(CASelectrons,MAXPTR,ptr1);
2760 nCASorb = str_nelem(CASorbitals,MAXPTR,ptr2);
2761 snew(ir->opts.CASelectrons,nr);
2762 snew(ir->opts.CASorbitals,nr);
2764 ir->opts.CASelectrons[i]= strtol(ptr1[i],NULL,10);
2765 ir->opts.CASorbitals[i] = strtol(ptr2[i],NULL,10);
2767 /* special optimization options */
2769 nbOPT = str_nelem(bOPT,MAXPTR,ptr1);
2770 nbTS = str_nelem(bTS,MAXPTR,ptr2);
2771 snew(ir->opts.bOPT,nr);
2772 snew(ir->opts.bTS,nr);
2774 ir->opts.bOPT[i] = (gmx_strncasecmp(ptr1[i],"Y",1)==0);
2775 ir->opts.bTS[i] = (gmx_strncasecmp(ptr2[i],"Y",1)==0);
2777 nSAon = str_nelem(SAon,MAXPTR,ptr1);
2778 nSAoff = str_nelem(SAoff,MAXPTR,ptr2);
2779 nSAsteps = str_nelem(SAsteps,MAXPTR,ptr3);
2780 snew(ir->opts.SAon,nr);
2781 snew(ir->opts.SAoff,nr);
2782 snew(ir->opts.SAsteps,nr);
2785 ir->opts.SAon[i] = strtod(ptr1[i],NULL);
2786 ir->opts.SAoff[i] = strtod(ptr2[i],NULL);
2787 ir->opts.SAsteps[i] = strtol(ptr3[i],NULL,10);
2789 /* end of QMMM input */
2792 for(i=0; (i<egcNR); i++) {
2793 fprintf(stderr,"%-16s has %d element(s):",gtypes[i],groups->grps[i].nr);
2794 for(j=0; (j<groups->grps[i].nr); j++)
2795 fprintf(stderr," %s",*(groups->grpname[groups->grps[i].nm_ind[j]]));
2796 fprintf(stderr,"\n");
2799 nr = groups->grps[egcENER].nr;
2800 snew(ir->opts.egp_flags,nr*nr);
2802 bExcl = do_egp_flag(ir,groups,"energygrp-excl",egpexcl,EGP_EXCL);
2803 if (bExcl && EEL_FULL(ir->coulombtype))
2804 warning(wi,"Can not exclude the lattice Coulomb energy between energy groups");
2806 bTable = do_egp_flag(ir,groups,"energygrp-table",egptable,EGP_TABLE);
2807 if (bTable && !(ir->vdwtype == evdwUSER) &&
2808 !(ir->coulombtype == eelUSER) && !(ir->coulombtype == eelPMEUSER) &&
2809 !(ir->coulombtype == eelPMEUSERSWITCH))
2810 gmx_fatal(FARGS,"Can only have energy group pair tables in combination with user tables for VdW and/or Coulomb");
2812 decode_cos(efield_x,&(ir->ex[XX]),FALSE);
2813 decode_cos(efield_xt,&(ir->et[XX]),TRUE);
2814 decode_cos(efield_y,&(ir->ex[YY]),FALSE);
2815 decode_cos(efield_yt,&(ir->et[YY]),TRUE);
2816 decode_cos(efield_z,&(ir->ex[ZZ]),FALSE);
2817 decode_cos(efield_zt,&(ir->et[ZZ]),TRUE);
2820 do_adress_index(ir->adress,groups,gnames,&(ir->opts),wi);
2822 for(i=0; (i<grps->nr); i++)
2832 static void check_disre(gmx_mtop_t *mtop)
2834 gmx_ffparams_t *ffparams;
2835 t_functype *functype;
2837 int i,ndouble,ftype;
2838 int label,old_label;
2840 if (gmx_mtop_ftype_count(mtop,F_DISRES) > 0) {
2841 ffparams = &mtop->ffparams;
2842 functype = ffparams->functype;
2843 ip = ffparams->iparams;
2846 for(i=0; i<ffparams->ntypes; i++) {
2847 ftype = functype[i];
2848 if (ftype == F_DISRES) {
2849 label = ip[i].disres.label;
2850 if (label == old_label) {
2851 fprintf(stderr,"Distance restraint index %d occurs twice\n",label);
2858 gmx_fatal(FARGS,"Found %d double distance restraint indices,\n"
2859 "probably the parameters for multiple pairs in one restraint "
2860 "are not identical\n",ndouble);
2864 static gmx_bool absolute_reference(t_inputrec *ir,gmx_mtop_t *sys,
2865 gmx_bool posres_only,
2869 gmx_mtop_ilistloop_t iloop;
2879 for(d=0; d<DIM; d++)
2881 AbsRef[d] = (d < ndof_com(ir) ? 0 : 1);
2883 /* Check for freeze groups */
2884 for(g=0; g<ir->opts.ngfrz; g++)
2886 for(d=0; d<DIM; d++)
2888 if (ir->opts.nFreeze[g][d] != 0)
2896 /* Check for position restraints */
2897 iloop = gmx_mtop_ilistloop_init(sys);
2898 while (gmx_mtop_ilistloop_next(iloop,&ilist,&nmol))
2901 (AbsRef[XX] == 0 || AbsRef[YY] == 0 || AbsRef[ZZ] == 0))
2903 for(i=0; i<ilist[F_POSRES].nr; i+=2)
2905 pr = &sys->ffparams.iparams[ilist[F_POSRES].iatoms[i]];
2906 for(d=0; d<DIM; d++)
2908 if (pr->posres.fcA[d] != 0)
2917 return (AbsRef[XX] != 0 && AbsRef[YY] != 0 && AbsRef[ZZ] != 0);
2920 void triple_check(const char *mdparin,t_inputrec *ir,gmx_mtop_t *sys,
2924 int i,m,g,nmol,npct;
2925 gmx_bool bCharge,bAcc;
2926 real gdt_max,*mgrp,mt;
2928 gmx_mtop_atomloop_block_t aloopb;
2929 gmx_mtop_atomloop_all_t aloop;
2932 char warn_buf[STRLEN];
2934 set_warning_line(wi,mdparin,-1);
2936 if (EI_DYNAMICS(ir->eI) && !EI_SD(ir->eI) && ir->eI != eiBD &&
2937 ir->comm_mode == ecmNO &&
2938 !(absolute_reference(ir,sys,FALSE,AbsRef) || ir->nsteps <= 10)) {
2939 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");
2942 /* Check for pressure coupling with absolute position restraints */
2943 if (ir->epc != epcNO && ir->refcoord_scaling == erscNO)
2945 absolute_reference(ir,sys,TRUE,AbsRef);
2947 for(m=0; m<DIM; m++)
2949 if (AbsRef[m] && norm2(ir->compress[m]) > 0)
2951 warning(wi,"You are using pressure coupling with absolute position restraints, this will give artifacts. Use the refcoord_scaling option.");
2959 aloopb = gmx_mtop_atomloop_block_init(sys);
2960 while (gmx_mtop_atomloop_block_next(aloopb,&atom,&nmol)) {
2961 if (atom->q != 0 || atom->qB != 0) {
2967 if (EEL_FULL(ir->coulombtype)) {
2969 "You are using full electrostatics treatment %s for a system without charges.\n"
2970 "This costs a lot of performance for just processing zeros, consider using %s instead.\n",
2971 EELTYPE(ir->coulombtype),EELTYPE(eelCUT));
2972 warning(wi,err_buf);
2975 if (ir->coulombtype == eelCUT && ir->rcoulomb > 0 && !ir->implicit_solvent) {
2977 "You are using a plain Coulomb cut-off, which might produce artifacts.\n"
2978 "You might want to consider using %s electrostatics.\n",
2980 warning_note(wi,err_buf);
2984 /* Generalized reaction field */
2985 if (ir->opts.ngtc == 0) {
2986 sprintf(err_buf,"No temperature coupling while using coulombtype %s",
2988 CHECK(ir->coulombtype == eelGRF);
2991 sprintf(err_buf,"When using coulombtype = %s"
2992 " ref-t for temperature coupling should be > 0",
2994 CHECK((ir->coulombtype == eelGRF) && (ir->opts.ref_t[0] <= 0));
2997 if (ir->eI == eiSD1 &&
2998 (gmx_mtop_ftype_count(sys,F_CONSTR) > 0 ||
2999 gmx_mtop_ftype_count(sys,F_SETTLE) > 0))
3001 sprintf(warn_buf,"With constraints integrator %s is less accurate, consider using %s instead",ei_names[ir->eI],ei_names[eiSD2]);
3002 warning_note(wi,warn_buf);
3006 for(i=0; (i<sys->groups.grps[egcACC].nr); i++) {
3007 for(m=0; (m<DIM); m++) {
3008 if (fabs(ir->opts.acc[i][m]) > 1e-6) {
3015 snew(mgrp,sys->groups.grps[egcACC].nr);
3016 aloop = gmx_mtop_atomloop_all_init(sys);
3017 while (gmx_mtop_atomloop_all_next(aloop,&i,&atom)) {
3018 mgrp[ggrpnr(&sys->groups,egcACC,i)] += atom->m;
3021 for(i=0; (i<sys->groups.grps[egcACC].nr); i++) {
3022 for(m=0; (m<DIM); m++)
3023 acc[m] += ir->opts.acc[i][m]*mgrp[i];
3026 for(m=0; (m<DIM); m++) {
3027 if (fabs(acc[m]) > 1e-6) {
3028 const char *dim[DIM] = { "X", "Y", "Z" };
3030 "Net Acceleration in %s direction, will %s be corrected\n",
3031 dim[m],ir->nstcomm != 0 ? "" : "not");
3032 if (ir->nstcomm != 0 && m < ndof_com(ir)) {
3034 for (i=0; (i<sys->groups.grps[egcACC].nr); i++)
3035 ir->opts.acc[i][m] -= acc[m];
3042 if (ir->efep != efepNO && ir->fepvals->sc_alpha != 0 &&
3043 !gmx_within_tol(sys->ffparams.reppow,12.0,10*GMX_DOUBLE_EPS)) {
3044 gmx_fatal(FARGS,"Soft-core interactions are only supported with VdW repulsion power 12");
3047 if (ir->ePull != epullNO) {
3048 if (ir->pull->grp[0].nat == 0) {
3049 absolute_reference(ir,sys,FALSE,AbsRef);
3050 for(m=0; m<DIM; m++) {
3051 if (ir->pull->dim[m] && !AbsRef[m]) {
3052 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.");
3058 if (ir->pull->eGeom == epullgDIRPBC) {
3059 for(i=0; i<3; i++) {
3060 for(m=0; m<=i; m++) {
3061 if ((ir->epc != epcNO && ir->compress[i][m] != 0) ||
3062 ir->deform[i][m] != 0) {
3063 for(g=1; g<ir->pull->ngrp; g++) {
3064 if (ir->pull->grp[g].vec[m] != 0) {
3065 gmx_fatal(FARGS,"Can not have dynamic box while using pull geometry '%s' (dim %c)",EPULLGEOM(ir->pull->eGeom),'x'+m);
3077 void double_check(t_inputrec *ir,matrix box,gmx_bool bConstr,warninp_t wi)
3081 char warn_buf[STRLEN];
3084 ptr = check_box(ir->ePBC,box);
3086 warning_error(wi,ptr);
3089 if (bConstr && ir->eConstrAlg == econtSHAKE) {
3090 if (ir->shake_tol <= 0.0) {
3091 sprintf(warn_buf,"ERROR: shake-tol must be > 0 instead of %g\n",
3093 warning_error(wi,warn_buf);
3096 if (IR_TWINRANGE(*ir) && ir->nstlist > 1) {
3097 sprintf(warn_buf,"With twin-range cut-off's and SHAKE the virial and the pressure are incorrect.");
3098 if (ir->epc == epcNO) {
3099 warning(wi,warn_buf);
3101 warning_error(wi,warn_buf);
3106 if( (ir->eConstrAlg == econtLINCS) && bConstr) {
3107 /* If we have Lincs constraints: */
3108 if(ir->eI==eiMD && ir->etc==etcNO &&
3109 ir->eConstrAlg==econtLINCS && ir->nLincsIter==1) {
3110 sprintf(warn_buf,"For energy conservation with LINCS, lincs_iter should be 2 or larger.\n");
3111 warning_note(wi,warn_buf);
3114 if ((ir->eI == eiCG || ir->eI == eiLBFGS) && (ir->nProjOrder<8)) {
3115 sprintf(warn_buf,"For accurate %s with LINCS constraints, lincs-order should be 8 or more.",ei_names[ir->eI]);
3116 warning_note(wi,warn_buf);
3118 if (ir->epc==epcMTTK) {
3119 warning_error(wi,"MTTK not compatible with lincs -- use shake instead.");
3123 if (ir->LincsWarnAngle > 90.0) {
3124 sprintf(warn_buf,"lincs-warnangle can not be larger than 90 degrees, setting it to 90.\n");
3125 warning(wi,warn_buf);
3126 ir->LincsWarnAngle = 90.0;
3129 if (ir->ePBC != epbcNONE) {
3130 if (ir->nstlist == 0) {
3131 warning(wi,"With nstlist=0 atoms are only put into the box at step 0, therefore drifting atoms might cause the simulation to crash.");
3133 bTWIN = (ir->rlistlong > ir->rlist);
3134 if (ir->ns_type == ensGRID) {
3135 if (sqr(ir->rlistlong) >= max_cutoff2(ir->ePBC,box)) {
3136 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",
3137 bTWIN ? (ir->rcoulomb==ir->rlistlong ? "rcoulomb" : "rvdw"):"rlist");
3138 warning_error(wi,warn_buf);
3141 min_size = min(box[XX][XX],min(box[YY][YY],box[ZZ][ZZ]));
3142 if (2*ir->rlistlong >= min_size) {
3143 sprintf(warn_buf,"ERROR: One of the box lengths is smaller than twice the cut-off length. Increase the box size or decrease rlist.");
3144 warning_error(wi,warn_buf);
3146 fprintf(stderr,"Grid search might allow larger cut-off's than simple search with triclinic boxes.");
3152 void check_chargegroup_radii(const gmx_mtop_t *mtop,const t_inputrec *ir,
3156 real rvdw1,rvdw2,rcoul1,rcoul2;
3157 char warn_buf[STRLEN];
3159 calc_chargegroup_radii(mtop,x,&rvdw1,&rvdw2,&rcoul1,&rcoul2);
3163 printf("Largest charge group radii for Van der Waals: %5.3f, %5.3f nm\n",
3168 printf("Largest charge group radii for Coulomb: %5.3f, %5.3f nm\n",
3174 if (rvdw1 + rvdw2 > ir->rlist ||
3175 rcoul1 + rcoul2 > ir->rlist)
3177 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);
3178 warning(wi,warn_buf);
3182 /* Here we do not use the zero at cut-off macro,
3183 * since user defined interactions might purposely
3184 * not be zero at the cut-off.
3186 if (EVDW_IS_ZERO_AT_CUTOFF(ir->vdwtype) &&
3187 rvdw1 + rvdw2 > ir->rlist - ir->rvdw)
3189 sprintf(warn_buf,"The sum of the two largest charge group radii (%f) is larger than rlist (%f) - rvdw (%f)\n",
3191 ir->rlist,ir->rvdw);
3194 warning(wi,warn_buf);
3198 warning_note(wi,warn_buf);
3201 if (EEL_IS_ZERO_AT_CUTOFF(ir->coulombtype) &&
3202 rcoul1 + rcoul2 > ir->rlistlong - ir->rcoulomb)
3204 sprintf(warn_buf,"The sum of the two largest charge group radii (%f) is larger than %s (%f) - rcoulomb (%f)\n",
3206 ir->rlistlong > ir->rlist ? "rlistlong" : "rlist",
3207 ir->rlistlong,ir->rcoulomb);
3210 warning(wi,warn_buf);
3214 warning_note(wi,warn_buf);