src/kernel/readir.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team,
   6  * check out http://www.gromacs.org for more information.
   7  * Copyright (c) 2012,2013, by the GROMACS development team, led by
   8  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   9  * others, as listed in the AUTHORS file in the top-level source
  10  * directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #ifdef HAVE_CONFIG_H
  39 #include <config.h>
  40 #endif
  41
  42 #include <ctype.h>
  43 #include <stdlib.h>
  44 #include <limits.h>
  45 #include "sysstuff.h"
  46 #include "smalloc.h"
  47 #include "typedefs.h"
  48 #include "physics.h"
  49 #include "names.h"
  50 #include "gmx_fatal.h"
  51 #include "macros.h"
  52 #include "index.h"
  53 #include "symtab.h"
  54 #include "string2.h"
  55 #include "readinp.h"
  56 #include "warninp.h"
  57 #include "readir.h"
  58 #include "toputil.h"
  59 #include "index.h"
  60 #include "network.h"
  61 #include "vec.h"
  62 #include "pbc.h"
  63 #include "mtop_util.h"
  64 #include "chargegroup.h"
  65 #include "inputrec.h"
  66
  67 #define MAXPTR 254
  68 #define NOGID  255
  69 #define MAXLAMBDAS 1024
  70
  71 /* Resource parameters
  72  * Do not change any of these until you read the instruction
  73  * in readinp.h. Some cpp's do not take spaces after the backslash
  74  * (like the c-shell), which will give you a very weird compiler
  75  * message.
  76  */
  77
  78 static char tcgrps[STRLEN],tau_t[STRLEN],ref_t[STRLEN],
  79   acc[STRLEN],accgrps[STRLEN],freeze[STRLEN],frdim[STRLEN],
  80   energy[STRLEN],user1[STRLEN],user2[STRLEN],vcm[STRLEN],xtc_grps[STRLEN],
  81   couple_moltype[STRLEN],orirefitgrp[STRLEN],egptable[STRLEN],egpexcl[STRLEN],
  82   wall_atomtype[STRLEN],wall_density[STRLEN],deform[STRLEN],QMMM[STRLEN];
  83 static char fep_lambda[efptNR][STRLEN];
  84 static char lambda_weights[STRLEN];
  85 static char **pull_grp;
  86 static char **rot_grp;
  87 static char anneal[STRLEN],anneal_npoints[STRLEN],
  88   anneal_time[STRLEN],anneal_temp[STRLEN];
  89 static char QMmethod[STRLEN],QMbasis[STRLEN],QMcharge[STRLEN],QMmult[STRLEN],
  90   bSH[STRLEN],CASorbitals[STRLEN], CASelectrons[STRLEN],SAon[STRLEN],
  91   SAoff[STRLEN],SAsteps[STRLEN],bTS[STRLEN],bOPT[STRLEN];
  92 static char efield_x[STRLEN],efield_xt[STRLEN],efield_y[STRLEN],
  93   efield_yt[STRLEN],efield_z[STRLEN],efield_zt[STRLEN];
  94
  95 enum {
  96     egrptpALL,         /* All particles have to be a member of a group.     */
  97     egrptpALL_GENREST, /* A rest group with name is generated for particles *
  98                         * that are not part of any group.                   */
  99     egrptpPART,        /* As egrptpALL_GENREST, but no name is generated    *
 100                         * for the rest group.                               */
 101     egrptpONE          /* Merge all selected groups into one group,         *
 102                         * make a rest group for the remaining particles.    */
 103 };
 104
 105
 106 void init_ir(t_inputrec *ir, t_gromppopts *opts)
 107 {
 108   snew(opts->include,STRLEN);
 109   snew(opts->define,STRLEN);
 110   snew(ir->fepvals,1);
 111   snew(ir->expandedvals,1);
 112   snew(ir->simtempvals,1);
 113 }
 114
 115 static void GetSimTemps(int ntemps, t_simtemp *simtemp, double *temperature_lambdas)
 116 {
 117
 118     int i;
 119
 120     for (i=0;i<ntemps;i++)
 121     {
 122         /* simple linear scaling -- allows more control */
 123         if (simtemp->eSimTempScale == esimtempLINEAR)
 124         {
 125             simtemp->temperatures[i] = simtemp->simtemp_low + (simtemp->simtemp_high-simtemp->simtemp_low)*temperature_lambdas[i];
 126         }
 127         else if (simtemp->eSimTempScale == esimtempGEOMETRIC)  /* should give roughly equal acceptance for constant heat capacity . . . */
 128         {
 129             simtemp->temperatures[i] = simtemp->simtemp_low * pow(simtemp->simtemp_high/simtemp->simtemp_low,(1.0*i)/(ntemps-1));
 130         }
 131         else if (simtemp->eSimTempScale == esimtempEXPONENTIAL)
 132         {
 133             simtemp->temperatures[i] = simtemp->simtemp_low + (simtemp->simtemp_high-simtemp->simtemp_low)*((exp(temperature_lambdas[i])-1)/(exp(1.0)-1));
 134         }
 135         else
 136         {
 137             char errorstr[128];
 138             sprintf(errorstr,"eSimTempScale=%d not defined",simtemp->eSimTempScale);
 139             gmx_fatal(FARGS,errorstr);
 140         }
 141     }
 142 }
 143
 144
 145
 146 static void _low_check(gmx_bool b,char *s,warninp_t wi)
 147 {
 148     if (b)
 149     {
 150         warning_error(wi,s);
 151     }
 152 }
 153
 154 static void check_nst(const char *desc_nst,int nst,
 155                       const char *desc_p,int *p,
 156                       warninp_t wi)
 157 {
 158     char buf[STRLEN];
 159
 160     if (*p > 0 && *p % nst != 0)
 161     {
 162         /* Round up to the next multiple of nst */
 163         *p = ((*p)/nst + 1)*nst;
 164         sprintf(buf,"%s should be a multiple of %s, changing %s to %d\n",
 165                 desc_p,desc_nst,desc_p,*p);
 166         warning(wi,buf);
 167     }
 168 }
 169
 170 static gmx_bool ir_NVE(const t_inputrec *ir)
 171 {
 172     return ((ir->eI == eiMD || EI_VV(ir->eI)) && ir->etc == etcNO);
 173 }
 174
 175 static int lcd(int n1,int n2)
 176 {
 177     int d,i;
 178
 179     d = 1;
 180     for(i=2; (i<=n1 && i<=n2); i++)
 181     {
 182         if (n1 % i == 0 && n2 % i == 0)
 183         {
 184             d = i;
 185         }
 186     }
 187
 188   return d;
 189 }
 190
 191 static void process_interaction_modifier(const t_inputrec *ir,int *eintmod)
 192 {
 193     if (*eintmod == eintmodPOTSHIFT_VERLET)
 194     {
 195         if (ir->cutoff_scheme == ecutsVERLET)
 196         {
 197             *eintmod = eintmodPOTSHIFT;
 198         }
 199         else
 200         {
 201             *eintmod = eintmodNONE;
 202         }
 203     }
 204 }
 205
 206 void check_ir(const char *mdparin,t_inputrec *ir, t_gromppopts *opts,
 207               warninp_t wi)
 208 /* Check internal consistency */
 209 {
 210     /* Strange macro: first one fills the err_buf, and then one can check
 211      * the condition, which will print the message and increase the error
 212      * counter.
 213      */
 214 #define CHECK(b) _low_check(b,err_buf,wi)
 215     char err_buf[256],warn_buf[STRLEN];
 216     int i,j;
 217     int  ns_type=0;
 218     real dt_coupl=0;
 219     real dt_pcoupl;
 220     int  nstcmin;
 221     t_lambda *fep = ir->fepvals;
 222     t_expanded *expand = ir->expandedvals;
 223
 224   set_warning_line(wi,mdparin,-1);
 225
 226     /* BASIC CUT-OFF STUFF */
 227     if (ir->rcoulomb < 0)
 228     {
 229         warning_error(wi,"rcoulomb should be >= 0");
 230     }
 231     if (ir->rvdw < 0)
 232     {
 233         warning_error(wi,"rvdw should be >= 0");
 234     }
 235     if (ir->rlist < 0 &&
 236         !(ir->cutoff_scheme == ecutsVERLET && ir->verletbuf_drift > 0))
 237     {
 238         warning_error(wi,"rlist should be >= 0");
 239     }
 240
 241     process_interaction_modifier(ir,&ir->coulomb_modifier);
 242     process_interaction_modifier(ir,&ir->vdw_modifier);
 243
 244     if (ir->cutoff_scheme == ecutsGROUP)
 245     {
 246         /* BASIC CUT-OFF STUFF */
 247         if (ir->rlist == 0 ||
 248             !((EEL_MIGHT_BE_ZERO_AT_CUTOFF(ir->coulombtype) && ir->rcoulomb > ir->rlist) ||
 249               (EVDW_MIGHT_BE_ZERO_AT_CUTOFF(ir->vdwtype)    && ir->rvdw     > ir->rlist))) {
 250             /* No switched potential and/or no twin-range:
 251              * we can set the long-range cut-off to the maximum of the other cut-offs.
 252              */
 253             ir->rlistlong = max_cutoff(ir->rlist,max_cutoff(ir->rvdw,ir->rcoulomb));
 254         }
 255         else if (ir->rlistlong < 0)
 256         {
 257             ir->rlistlong = max_cutoff(ir->rlist,max_cutoff(ir->rvdw,ir->rcoulomb));
 258             sprintf(warn_buf,"rlistlong was not set, setting it to %g (no buffer)",
 259                     ir->rlistlong);
 260             warning(wi,warn_buf);
 261         }
 262         if (ir->rlistlong == 0 && ir->ePBC != epbcNONE)
 263         {
 264             warning_error(wi,"Can not have an infinite cut-off with PBC");
 265         }
 266         if (ir->rlistlong > 0 && (ir->rlist == 0 || ir->rlistlong < ir->rlist))
 267         {
 268             warning_error(wi,"rlistlong can not be shorter than rlist");
 269         }
 270         if (IR_TWINRANGE(*ir) && ir->nstlist <= 0)
 271         {
 272             warning_error(wi,"Can not have nstlist<=0 with twin-range interactions");
 273         }
 274     }
 275
 276     if(ir->rlistlong == ir->rlist)
 277     {
 278         ir->nstcalclr = 0;
 279     }
 280     else if(ir->rlistlong>ir->rlist && ir->nstcalclr==0)
 281     {
 282         warning_error(wi,"With different cutoffs for electrostatics and VdW, nstcalclr must be -1 or a positive number");
 283     }
 284
 285     if (ir->cutoff_scheme == ecutsVERLET)
 286     {
 287         real rc_max;
 288
 289         /* Normal Verlet type neighbor-list, currently only limited feature support */
 290         if (inputrec2nboundeddim(ir) < 3)
 291         {
 292             warning_error(wi,"With Verlet lists only full pbc or pbc=xy with walls is supported");
 293         }
 294         if (ir->rcoulomb != ir->rvdw)
 295         {
 296             warning_error(wi,"With Verlet lists rcoulomb!=rvdw is not supported");
 297         }
 298         if (ir->vdwtype != evdwCUT)
 299         {
 300             warning_error(wi,"With Verlet lists only cut-off LJ interactions are supported");
 301         }
 302         if (!(ir->coulombtype == eelCUT ||
 303               (EEL_RF(ir->coulombtype) && ir->coulombtype != eelRF_NEC) ||
 304               EEL_PME(ir->coulombtype) || ir->coulombtype == eelEWALD))
 305         {
 306             warning_error(wi,"With Verlet lists only cut-off, reaction-field, PME and Ewald electrostatics are supported");
 307         }
 308
 309         if (ir->nstlist <= 0)
 310         {
 311              warning_error(wi,"With Verlet lists nstlist should be larger than 0");
 312         }
 313
 314         if (ir->nstlist < 10)
 315         {
 316             warning_note(wi,"With Verlet lists the optimal nstlist is >= 10, with GPUs >= 20. Note that with the Verlet scheme, nstlist has no effect on the accuracy of your simulation.");
 317         }
 318
 319         rc_max = max(ir->rvdw,ir->rcoulomb);
 320
 321         if (ir->verletbuf_drift <= 0)
 322         {
 323             if (ir->verletbuf_drift == 0)
 324             {
 325                 warning_error(wi,"Can not have an energy drift of exactly 0");
 326             }
 327
 328             if (ir->rlist < rc_max)
 329             {
 330                 warning_error(wi,"With verlet lists rlist can not be smaller than rvdw or rcoulomb");
 331             }
 332
 333             if (ir->rlist == rc_max && ir->nstlist > 1)
 334             {
 335                 warning_note(wi,"rlist is equal to rvdw and/or rcoulomb: there is no explicit Verlet buffer. The cluster pair list does have a buffering effect, but choosing a larger rlist might be necessary for good energy conservation.");
 336             }
 337         }
 338         else
 339         {
 340             if (ir->rlist > rc_max)
 341             {
 342                 warning_note(wi,"You have set rlist larger than the interaction cut-off, but you also have verlet-buffer-drift > 0. Will set rlist using verlet-buffer-drift.");
 343             }
 344
 345             if (ir->nstlist == 1)
 346             {
 347                 /* No buffer required */
 348                 ir->rlist = rc_max;
 349             }
 350             else
 351             {
 352                 if (EI_DYNAMICS(ir->eI))
 353                 {
 354                     if (EI_MD(ir->eI) && ir->etc == etcNO)
 355                     {
 356                         warning_error(wi,"Temperature coupling is required for calculating rlist using the energy drift with verlet-buffer-drift > 0. Either use temperature coupling or set rlist yourself together with verlet-buffer-drift = -1.");
 357                     }
 358
 359                     if (inputrec2nboundeddim(ir) < 3)
 360                     {
 361                         warning_error(wi,"The box volume is required for calculating rlist from the energy drift with verlet-buffer-drift > 0. You are using at least one unbounded dimension, so no volume can be computed. Either use a finite box, or set rlist yourself together with verlet-buffer-drift = -1.");
 362                     }
 363                     /* Set rlist temporarily so we can continue processing */
 364                     ir->rlist = rc_max;
 365                 }
 366                 else
 367                 {
 368                     /* Set the buffer to 5% of the cut-off */
 369                     ir->rlist = 1.05*rc_max;
 370                 }
 371             }
 372         }
 373
 374         /* No twin-range calculations with Verlet lists */
 375         ir->rlistlong = ir->rlist;
 376     }
 377
 378     if(ir->nstcalclr==-1)
 379     {
 380         /* if rlist=rlistlong, this will later be changed to nstcalclr=0 */
 381         ir->nstcalclr = ir->nstlist;
 382     }
 383     else if(ir->nstcalclr>0)
 384     {
 385         if(ir->nstlist>0 && (ir->nstlist % ir->nstcalclr != 0))
 386         {
 387             warning_error(wi,"nstlist must be evenly divisible by nstcalclr. Use nstcalclr = -1 to automatically follow nstlist");
 388         }
 389     }
 390     else if(ir->nstcalclr<-1)
 391     {
 392         warning_error(wi,"nstcalclr must be a positive number (divisor of nstcalclr), or -1 to follow nstlist.");
 393     }
 394
 395     if(EEL_PME(ir->coulombtype) && ir->rcoulomb > ir->rvdw && ir->nstcalclr>1)
 396     {
 397         warning_error(wi,"When used with PME, the long-range component of twin-range interactions must be updated every step (nstcalclr)");
 398     }
 399
 400     /* GENERAL INTEGRATOR STUFF */
 401     if (!(ir->eI == eiMD || EI_VV(ir->eI)))
 402     {
 403         ir->etc = etcNO;
 404     }
 405     if (ir->eI == eiVVAK) {
 406         sprintf(warn_buf,"Integrator method %s is implemented primarily for validation purposes; for molecular dynamics, you should probably be using %s or %s",ei_names[eiVVAK],ei_names[eiMD],ei_names[eiVV]);
 407         warning_note(wi,warn_buf);
 408     }
 409     if (!EI_DYNAMICS(ir->eI))
 410     {
 411         ir->epc = epcNO;
 412     }
 413     if (EI_DYNAMICS(ir->eI))
 414     {
 415         if (ir->nstcalcenergy < 0)
 416         {
 417             ir->nstcalcenergy = ir_optimal_nstcalcenergy(ir);
 418             if (ir->nstenergy != 0 && ir->nstenergy < ir->nstcalcenergy)
 419             {
 420                 /* nstcalcenergy larger than nstener does not make sense.
 421                  * We ideally want nstcalcenergy=nstener.
 422                  */
 423                 if (ir->nstlist > 0)
 424                 {
 425                     ir->nstcalcenergy = lcd(ir->nstenergy,ir->nstlist);
 426                 }
 427                 else
 428                 {
 429                     ir->nstcalcenergy = ir->nstenergy;
 430                 }
 431             }
 432         }
 433         else if ( (ir->nstenergy > 0 && ir->nstcalcenergy > ir->nstenergy) ||
 434                   (ir->efep != efepNO && ir->fepvals->nstdhdl > 0 &&
 435                    (ir->nstcalcenergy > ir->fepvals->nstdhdl) ) )
 436
 437         {
 438             const char *nsten="nstenergy";
 439             const char *nstdh="nstdhdl";
 440             const char *min_name=nsten;
 441             int min_nst=ir->nstenergy;
 442
 443             /* find the smallest of ( nstenergy, nstdhdl ) */
 444             if (ir->efep != efepNO && ir->fepvals->nstdhdl > 0 &&
 445                 (ir->fepvals->nstdhdl < ir->nstenergy) )
 446             {
 447                 min_nst=ir->fepvals->nstdhdl;
 448                 min_name=nstdh;
 449             }
 450             /* If the user sets nstenergy small, we should respect that */
 451             sprintf(warn_buf,
 452                     "Setting nstcalcenergy (%d) equal to %s (%d)",
 453                     ir->nstcalcenergy,min_name, min_nst);
 454             warning_note(wi,warn_buf);
 455             ir->nstcalcenergy = min_nst;
 456         }
 457
 458         if (ir->epc != epcNO)
 459         {
 460             if (ir->nstpcouple < 0)
 461             {
 462                 ir->nstpcouple = ir_optimal_nstpcouple(ir);
 463             }
 464         }
 465         if (IR_TWINRANGE(*ir))
 466         {
 467             check_nst("nstlist",ir->nstlist,
 468                       "nstcalcenergy",&ir->nstcalcenergy,wi);
 469             if (ir->epc != epcNO)
 470             {
 471                 check_nst("nstlist",ir->nstlist,
 472                           "nstpcouple",&ir->nstpcouple,wi);
 473             }
 474         }
 475
 476         if (ir->nstcalcenergy > 0)
 477         {
 478             if (ir->efep != efepNO)
 479             {
 480                 /* nstdhdl should be a multiple of nstcalcenergy */
 481                 check_nst("nstcalcenergy",ir->nstcalcenergy,
 482                           "nstdhdl",&ir->fepvals->nstdhdl,wi);
 483                 /* nstexpanded should be a multiple of nstcalcenergy */
 484                 check_nst("nstcalcenergy",ir->nstcalcenergy,
 485                           "nstexpanded",&ir->expandedvals->nstexpanded,wi);
 486             }
 487             /* for storing exact averages nstenergy should be
 488              * a multiple of nstcalcenergy
 489              */
 490             check_nst("nstcalcenergy",ir->nstcalcenergy,
 491                       "nstenergy",&ir->nstenergy,wi);
 492         }
 493     }
 494
 495   /* LD STUFF */
 496   if ((EI_SD(ir->eI) || ir->eI == eiBD) &&
 497       ir->bContinuation && ir->ld_seed != -1) {
 498       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)");
 499   }
 500
 501   /* TPI STUFF */
 502   if (EI_TPI(ir->eI)) {
 503     sprintf(err_buf,"TPI only works with pbc = %s",epbc_names[epbcXYZ]);
 504     CHECK(ir->ePBC != epbcXYZ);
 505     sprintf(err_buf,"TPI only works with ns = %s",ens_names[ensGRID]);
 506     CHECK(ir->ns_type != ensGRID);
 507     sprintf(err_buf,"with TPI nstlist should be larger than zero");
 508     CHECK(ir->nstlist <= 0);
 509     sprintf(err_buf,"TPI does not work with full electrostatics other than PME");
 510     CHECK(EEL_FULL(ir->coulombtype) && !EEL_PME(ir->coulombtype));
 511   }
 512
 513   /* SHAKE / LINCS */
 514   if ( (opts->nshake > 0) && (opts->bMorse) ) {
 515       sprintf(warn_buf,
 516               "Using morse bond-potentials while constraining bonds is useless");
 517       warning(wi,warn_buf);
 518   }
 519
 520   if ((EI_SD(ir->eI) || ir->eI == eiBD) &&
 521       ir->bContinuation && ir->ld_seed != -1) {
 522       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)");
 523   }
 524   /* verify simulated tempering options */
 525
 526   if (ir->bSimTemp) {
 527       gmx_bool bAllTempZero = TRUE;
 528       for (i=0;i<fep->n_lambda;i++)
 529       {
 530           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]);
 531           CHECK((fep->all_lambda[efptTEMPERATURE][i] < 0) || (fep->all_lambda[efptTEMPERATURE][i] > 1));
 532           if (fep->all_lambda[efptTEMPERATURE][i] > 0)
 533           {
 534               bAllTempZero = FALSE;
 535           }
 536       }
 537       sprintf(err_buf,"if simulated tempering is on, temperature-lambdas may not be all zero");
 538       CHECK(bAllTempZero==TRUE);
 539
 540       sprintf(err_buf,"Simulated tempering is currently only compatible with md-vv");
 541       CHECK(ir->eI != eiVV);
 542
 543       /* check compatability of the temperature coupling with simulated tempering */
 544
 545       if (ir->etc == etcNOSEHOOVER) {
 546           sprintf(warn_buf,"Nose-Hoover based temperature control such as [%s] my not be entirelyconsistent with simulated tempering",etcoupl_names[ir->etc]);
 547           warning_note(wi,warn_buf);
 548       }
 549
 550       /* check that the temperatures make sense */
 551
 552       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);
 553       CHECK(ir->simtempvals->simtemp_high <= ir->simtempvals->simtemp_low);
 554
 555       sprintf(err_buf,"Higher simulated tempering temperature (%g) must be >= zero",ir->simtempvals->simtemp_high);
 556       CHECK(ir->simtempvals->simtemp_high <= 0);
 557
 558       sprintf(err_buf,"Lower simulated tempering temperature (%g) must be >= zero",ir->simtempvals->simtemp_low);
 559       CHECK(ir->simtempvals->simtemp_low <= 0);
 560   }
 561
 562   /* verify free energy options */
 563
 564   if (ir->efep != efepNO) {
 565       fep = ir->fepvals;
 566       sprintf(err_buf,"The soft-core power is %d and can only be 1 or 2",
 567               fep->sc_power);
 568       CHECK(fep->sc_alpha!=0 && fep->sc_power!=1 && fep->sc_power!=2);
 569
 570       sprintf(err_buf,"The soft-core sc-r-power is %d and can only be 6 or 48",
 571               (int)fep->sc_r_power);
 572       CHECK(fep->sc_alpha!=0 && fep->sc_r_power!=6.0 && fep->sc_r_power!=48.0);
 573
 574       /* check validity of options */
 575       if (fep->n_lambda > 0 && ir->rlist < max(ir->rvdw,ir->rcoulomb))
 576       {
 577           sprintf(warn_buf,
 578                   "For foreign lambda free energy differences it is assumed that the soft-core interactions have no effect beyond the neighborlist cut-off");
 579           warning(wi,warn_buf);
 580       }
 581
 582       sprintf(err_buf,"Can't use postive delta-lambda (%g) if initial state/lambda does not start at zero",fep->delta_lambda);
 583       CHECK(fep->delta_lambda > 0 && ((fep->init_fep_state > 0) ||  (fep->init_lambda > 0)));
 584
 585       sprintf(err_buf,"Can't use postive delta-lambda (%g) with expanded ensemble simulations",fep->delta_lambda);
 586       CHECK(fep->delta_lambda > 0 && (ir->efep == efepEXPANDED));
 587
 588       sprintf(err_buf,"Free-energy not implemented for Ewald");
 589       CHECK(ir->coulombtype==eelEWALD);
 590
 591       /* check validty of lambda inputs */
 592       if (fep->n_lambda == 0)
 593       {
 594           /* Clear output in case of no states:*/
 595           sprintf(err_buf,"init-lambda-state set to %d: no lambda states are defined.",fep->init_fep_state);
 596           CHECK((fep->init_fep_state>=0) && (fep->n_lambda==0));
 597       }
 598       else
 599       {
 600           sprintf(err_buf,"initial thermodynamic state %d does not exist, only goes to %d",fep->init_fep_state,fep->n_lambda-1);
 601           CHECK((fep->init_fep_state >= fep->n_lambda));
 602       }
 603
 604       sprintf(err_buf,"Lambda state must be set, either with init-lambda-state or with init-lambda");
 605       CHECK((fep->init_fep_state < 0) && (fep->init_lambda <0));
 606
 607       sprintf(err_buf,"init-lambda=%g while init-lambda-state=%d. Lambda state must be set either with init-lambda-state or with init-lambda, but not both",
 608               fep->init_lambda, fep->init_fep_state);
 609       CHECK((fep->init_fep_state >= 0) && (fep->init_lambda >= 0));
 610
 611
 612
 613       if((fep->init_lambda >= 0) && (fep->delta_lambda == 0))
 614       {
 615           int n_lambda_terms;
 616           n_lambda_terms=0;
 617           for (i=0;i<efptNR;i++)
 618           {
 619               if (fep->separate_dvdl[i])
 620               {
 621                   n_lambda_terms++;
 622               }
 623           }
 624           if (n_lambda_terms > 1)
 625           {
 626               sprintf(warn_buf,"If lambda vector states (fep-lambdas, coul-lambdas etc.) are set, don't use init-lambda to set lambda state (except for slow growth). Use init-lambda-state instead.");
 627               warning(wi, warn_buf);
 628           }
 629
 630           if (n_lambda_terms < 2 && fep->n_lambda > 0)
 631           {
 632               warning_note(wi,
 633                            "init-lambda is deprecated for setting lambda state (except for slow growth). Use init-lambda-state instead.");
 634           }
 635       }
 636
 637       for (j=0;j<efptNR;j++)
 638       {
 639           for (i=0;i<fep->n_lambda;i++)
 640           {
 641               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]);
 642               CHECK((fep->all_lambda[j][i] < 0) || (fep->all_lambda[j][i] > 1));
 643           }
 644       }
 645
 646       if ((fep->sc_alpha>0) && (!fep->bScCoul))
 647       {
 648           for (i=0;i<fep->n_lambda;i++)
 649           {
 650               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],
 651                       fep->all_lambda[efptCOUL][i]);
 652               CHECK((fep->sc_alpha>0) &&
 653                     (((fep->all_lambda[efptCOUL][i] > 0.0) &&
 654                       (fep->all_lambda[efptCOUL][i] < 1.0)) &&
 655                      ((fep->all_lambda[efptVDW][i] > 0.0) &&
 656                       (fep->all_lambda[efptVDW][i] < 1.0))));
 657           }
 658       }
 659
 660       if ((fep->bScCoul) && (EEL_PME(ir->coulombtype)))
 661       {
 662           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.");
 663           warning(wi, warn_buf);
 664       }
 665
 666       /*  Free Energy Checks -- In an ideal world, slow growth and FEP would
 667           be treated differently, but that's the next step */
 668
 669       for (i=0;i<efptNR;i++) {
 670           for (j=0;j<fep->n_lambda;j++) {
 671               sprintf(err_buf,"%s[%d] must be between 0 and 1",efpt_names[i],j);
 672               CHECK((fep->all_lambda[i][j] < 0) || (fep->all_lambda[i][j] > 1));
 673           }
 674       }
 675   }
 676
 677   if ((ir->bSimTemp) || (ir->efep == efepEXPANDED)) {
 678       fep = ir->fepvals;
 679       expand = ir->expandedvals;
 680
 681       /* checking equilibration of weights inputs for validity */
 682
 683       sprintf(err_buf,"weight-equil-number-all-lambda (%d) is ignored if lmc-weights-equil is not equal to %s",
 684               expand->equil_n_at_lam,elmceq_names[elmceqNUMATLAM]);
 685       CHECK((expand->equil_n_at_lam>0) && (expand->elmceq!=elmceqNUMATLAM));
 686
 687       sprintf(err_buf,"weight-equil-number-samples (%d) is ignored if lmc-weights-equil is not equal to %s",
 688               expand->equil_samples,elmceq_names[elmceqSAMPLES]);
 689       CHECK((expand->equil_samples>0) && (expand->elmceq!=elmceqSAMPLES));
 690
 691       sprintf(err_buf,"weight-equil-number-steps (%d) is ignored if lmc-weights-equil is not equal to %s",
 692               expand->equil_steps,elmceq_names[elmceqSTEPS]);
 693       CHECK((expand->equil_steps>0) && (expand->elmceq!=elmceqSTEPS));
 694
 695       sprintf(err_buf,"weight-equil-wl-delta (%d) is ignored if lmc-weights-equil is not equal to %s",
 696               expand->equil_samples,elmceq_names[elmceqWLDELTA]);
 697       CHECK((expand->equil_wl_delta>0) && (expand->elmceq!=elmceqWLDELTA));
 698
 699       sprintf(err_buf,"weight-equil-count-ratio (%f) is ignored if lmc-weights-equil is not equal to %s",
 700               expand->equil_ratio,elmceq_names[elmceqRATIO]);
 701       CHECK((expand->equil_ratio>0) && (expand->elmceq!=elmceqRATIO));
 702
 703       sprintf(err_buf,"weight-equil-number-all-lambda (%d) must be a positive integer if lmc-weights-equil=%s",
 704               expand->equil_n_at_lam,elmceq_names[elmceqNUMATLAM]);
 705       CHECK((expand->equil_n_at_lam<=0) && (expand->elmceq==elmceqNUMATLAM));
 706
 707       sprintf(err_buf,"weight-equil-number-samples (%d) must be a positive integer if lmc-weights-equil=%s",
 708               expand->equil_samples,elmceq_names[elmceqSAMPLES]);
 709       CHECK((expand->equil_samples<=0) && (expand->elmceq==elmceqSAMPLES));
 710
 711       sprintf(err_buf,"weight-equil-number-steps (%d) must be a positive integer if lmc-weights-equil=%s",
 712               expand->equil_steps,elmceq_names[elmceqSTEPS]);
 713       CHECK((expand->equil_steps<=0) && (expand->elmceq==elmceqSTEPS));
 714
 715       sprintf(err_buf,"weight-equil-wl-delta (%f) must be > 0 if lmc-weights-equil=%s",
 716               expand->equil_wl_delta,elmceq_names[elmceqWLDELTA]);
 717       CHECK((expand->equil_wl_delta<=0) && (expand->elmceq==elmceqWLDELTA));
 718
 719       sprintf(err_buf,"weight-equil-count-ratio (%f) must be > 0 if lmc-weights-equil=%s",
 720               expand->equil_ratio,elmceq_names[elmceqRATIO]);
 721       CHECK((expand->equil_ratio<=0) && (expand->elmceq==elmceqRATIO));
 722
 723       sprintf(err_buf,"lmc-weights-equil=%s only possible when lmc-stats = %s or lmc-stats %s",
 724               elmceq_names[elmceqWLDELTA],elamstats_names[elamstatsWL],elamstats_names[elamstatsWWL]);
 725       CHECK((expand->elmceq==elmceqWLDELTA) && (!EWL(expand->elamstats)));
 726
 727       sprintf(err_buf,"lmc-repeats (%d) must be greater than 0",expand->lmc_repeats);
 728       CHECK((expand->lmc_repeats <= 0));
 729       sprintf(err_buf,"minimum-var-min (%d) must be greater than 0",expand->minvarmin);
 730       CHECK((expand->minvarmin <= 0));
 731       sprintf(err_buf,"weight-c-range (%d) must be greater or equal to 0",expand->c_range);
 732       CHECK((expand->c_range < 0));
 733       sprintf(err_buf,"init-lambda-state (%d) must be zero if lmc-forced-nstart (%d)> 0 and lmc-move != 'no'",
 734               fep->init_fep_state, expand->lmc_forced_nstart);
 735       CHECK((fep->init_fep_state!=0) && (expand->lmc_forced_nstart>0) && (expand->elmcmove!=elmcmoveNO));
 736       sprintf(err_buf,"lmc-forced-nstart (%d) must not be negative",expand->lmc_forced_nstart);
 737       CHECK((expand->lmc_forced_nstart < 0));
 738       sprintf(err_buf,"init-lambda-state (%d) must be in the interval [0,number of lambdas)",fep->init_fep_state);
 739       CHECK((fep->init_fep_state < 0) || (fep->init_fep_state >= fep->n_lambda));
 740
 741       sprintf(err_buf,"init-wl-delta (%f) must be greater than or equal to 0",expand->init_wl_delta);
 742       CHECK((expand->init_wl_delta < 0));
 743       sprintf(err_buf,"wl-ratio (%f) must be between 0 and 1",expand->wl_ratio);
 744       CHECK((expand->wl_ratio <= 0) || (expand->wl_ratio >= 1));
 745       sprintf(err_buf,"wl-scale (%f) must be between 0 and 1",expand->wl_scale);
 746       CHECK((expand->wl_scale <= 0) || (expand->wl_scale >= 1));
 747
 748       /* if there is no temperature control, we need to specify an MC temperature */
 749       sprintf(err_buf,"If there is no temperature control, and lmc-mcmove!= 'no',mc_temperature must be set to a positive number");
 750       if (expand->nstTij > 0)
 751       {
 752           sprintf(err_buf,"nst-transition-matrix (%d) must be an integer multiple of nstlog (%d)",
 753                   expand->nstTij,ir->nstlog);
 754           CHECK((mod(expand->nstTij,ir->nstlog)!=0));
 755       }
 756   }
 757
 758   /* PBC/WALLS */
 759   sprintf(err_buf,"walls only work with pbc=%s",epbc_names[epbcXY]);
 760   CHECK(ir->nwall && ir->ePBC!=epbcXY);
 761
 762   /* VACUUM STUFF */
 763   if (ir->ePBC != epbcXYZ && ir->nwall != 2) {
 764     if (ir->ePBC == epbcNONE) {
 765       if (ir->epc != epcNO) {
 766           warning(wi,"Turning off pressure coupling for vacuum system");
 767           ir->epc = epcNO;
 768       }
 769     } else {
 770       sprintf(err_buf,"Can not have pressure coupling with pbc=%s",
 771               epbc_names[ir->ePBC]);
 772       CHECK(ir->epc != epcNO);
 773     }
 774     sprintf(err_buf,"Can not have Ewald with pbc=%s",epbc_names[ir->ePBC]);
 775     CHECK(EEL_FULL(ir->coulombtype));
 776
 777     sprintf(err_buf,"Can not have dispersion correction with pbc=%s",
 778             epbc_names[ir->ePBC]);
 779     CHECK(ir->eDispCorr != edispcNO);
 780   }
 781
 782   if (ir->rlist == 0.0) {
 783     sprintf(err_buf,"can only have neighborlist cut-off zero (=infinite)\n"
 784             "with coulombtype = %s or coulombtype = %s\n"
 785             "without periodic boundary conditions (pbc = %s) and\n"
 786             "rcoulomb and rvdw set to zero",
 787             eel_names[eelCUT],eel_names[eelUSER],epbc_names[epbcNONE]);
 788     CHECK(((ir->coulombtype != eelCUT) && (ir->coulombtype != eelUSER)) ||
 789           (ir->ePBC     != epbcNONE) ||
 790           (ir->rcoulomb != 0.0)      || (ir->rvdw != 0.0));
 791
 792     if (ir->nstlist < 0) {
 793         warning_error(wi,"Can not have heuristic neighborlist updates without cut-off");
 794     }
 795     if (ir->nstlist > 0) {
 796         warning_note(wi,"Simulating without cut-offs is usually (slightly) faster with nstlist=0, nstype=simple and particle decomposition");
 797     }
 798   }
 799
 800   /* COMM STUFF */
 801   if (ir->nstcomm == 0) {
 802     ir->comm_mode = ecmNO;
 803   }
 804   if (ir->comm_mode != ecmNO) {
 805     if (ir->nstcomm < 0) {
 806         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");
 807       ir->nstcomm = abs(ir->nstcomm);
 808     }
 809
 810     if (ir->nstcalcenergy > 0 && ir->nstcomm < ir->nstcalcenergy) {
 811         warning_note(wi,"nstcomm < nstcalcenergy defeats the purpose of nstcalcenergy, setting nstcomm to nstcalcenergy");
 812         ir->nstcomm = ir->nstcalcenergy;
 813     }
 814
 815     if (ir->comm_mode == ecmANGULAR) {
 816       sprintf(err_buf,"Can not remove the rotation around the center of mass with periodic molecules");
 817       CHECK(ir->bPeriodicMols);
 818       if (ir->ePBC != epbcNONE)
 819           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).");
 820     }
 821   }
 822
 823   if (EI_STATE_VELOCITY(ir->eI) && ir->ePBC == epbcNONE && ir->comm_mode != ecmANGULAR) {
 824       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.");
 825   }
 826
 827   sprintf(err_buf,"Twin-range neighbour searching (NS) with simple NS"
 828           " algorithm not implemented");
 829   CHECK(((ir->rcoulomb > ir->rlist) || (ir->rvdw > ir->rlist))
 830         && (ir->ns_type == ensSIMPLE));
 831
 832   /* TEMPERATURE COUPLING */
 833   if (ir->etc == etcYES)
 834     {
 835         ir->etc = etcBERENDSEN;
 836         warning_note(wi,"Old option for temperature coupling given: "
 837                      "changing \"yes\" to \"Berendsen\"\n");
 838     }
 839
 840     if ((ir->etc == etcNOSEHOOVER) || (ir->epc == epcMTTK))
 841     {
 842         if (ir->opts.nhchainlength < 1)
 843         {
 844             sprintf(warn_buf,"number of Nose-Hoover chains (currently %d) cannot be less than 1,reset to 1\n",ir->opts.nhchainlength);
 845             ir->opts.nhchainlength =1;
 846             warning(wi,warn_buf);
 847         }
 848
 849         if (ir->etc==etcNOSEHOOVER && !EI_VV(ir->eI) && ir->opts.nhchainlength > 1)
 850         {
 851             warning_note(wi,"leapfrog does not yet support Nose-Hoover chains, nhchainlength reset to 1");
 852             ir->opts.nhchainlength = 1;
 853         }
 854     }
 855     else
 856     {
 857         ir->opts.nhchainlength = 0;
 858     }
 859
 860     if (ir->eI == eiVVAK) {
 861         sprintf(err_buf,"%s implemented primarily for validation, and requires nsttcouple = 1 and nstpcouple = 1.",
 862                 ei_names[eiVVAK]);
 863         CHECK((ir->nsttcouple != 1) || (ir->nstpcouple != 1));
 864     }
 865
 866     if (ETC_ANDERSEN(ir->etc))
 867     {
 868         sprintf(err_buf,"%s temperature control not supported for integrator %s.",etcoupl_names[ir->etc],ei_names[ir->eI]);
 869         CHECK(!(EI_VV(ir->eI)));
 870
 871         for (i=0;i<ir->opts.ngtc;i++)
 872         {
 873             sprintf(err_buf,"all tau_t must currently be equal using Andersen temperature control, violated for group %d",i);
 874             CHECK(ir->opts.tau_t[0] != ir->opts.tau_t[i]);
 875             sprintf(err_buf,"all tau_t must be postive using Andersen temperature control, tau_t[%d]=%10.6f",
 876                     i,ir->opts.tau_t[i]);
 877             CHECK(ir->opts.tau_t[i]<0);
 878         }
 879         if (ir->nstcomm > 0 && (ir->etc == etcANDERSEN)) {
 880             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]);
 881             warning_note(wi,warn_buf);
 882         }
 883
 884         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]);
 885         CHECK(ir->nstcomm > 1 && (ir->etc == etcANDERSEN));
 886
 887         for (i=0;i<ir->opts.ngtc;i++)
 888         {
 889             int nsteps = (int)(ir->opts.tau_t[i]/ir->delta_t);
 890             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);
 891             CHECK((nsteps % ir->nstcomm) && (ir->etc == etcANDERSENMASSIVE));
 892         }
 893     }
 894     if (ir->etc == etcBERENDSEN)
 895     {
 896         sprintf(warn_buf,"The %s thermostat does not generate the correct kinetic energy distribution. You might want to consider using the %s thermostat.",
 897                 ETCOUPLTYPE(ir->etc),ETCOUPLTYPE(etcVRESCALE));
 898         warning_note(wi,warn_buf);
 899     }
 900
 901     if ((ir->etc==etcNOSEHOOVER || ETC_ANDERSEN(ir->etc))
 902         && ir->epc==epcBERENDSEN)
 903     {
 904         sprintf(warn_buf,"Using Berendsen pressure coupling invalidates the "
 905                 "true ensemble for the thermostat");
 906         warning(wi,warn_buf);
 907     }
 908
 909     /* PRESSURE COUPLING */
 910     if (ir->epc == epcISOTROPIC)
 911     {
 912         ir->epc = epcBERENDSEN;
 913         warning_note(wi,"Old option for pressure coupling given: "
 914                      "changing \"Isotropic\" to \"Berendsen\"\n");
 915     }
 916
 917     if (ir->epc != epcNO)
 918     {
 919         dt_pcoupl = ir->nstpcouple*ir->delta_t;
 920
 921         sprintf(err_buf,"tau-p must be > 0 instead of %g\n",ir->tau_p);
 922         CHECK(ir->tau_p <= 0);
 923
 924         if (ir->tau_p/dt_pcoupl < pcouple_min_integration_steps(ir->epc))
 925         {
 926             sprintf(warn_buf,"For proper integration of the %s barostat, tau-p (%g) should be at least %d times larger than nstpcouple*dt (%g)",
 927                     EPCOUPLTYPE(ir->epc),ir->tau_p,pcouple_min_integration_steps(ir->epc),dt_pcoupl);
 928             warning(wi,warn_buf);
 929         }
 930
 931         sprintf(err_buf,"compressibility must be > 0 when using pressure"
 932                 " coupling %s\n",EPCOUPLTYPE(ir->epc));
 933         CHECK(ir->compress[XX][XX] < 0 || ir->compress[YY][YY] < 0 ||
 934               ir->compress[ZZ][ZZ] < 0 ||
 935               (trace(ir->compress) == 0 && ir->compress[YY][XX] <= 0 &&
 936                ir->compress[ZZ][XX] <= 0 && ir->compress[ZZ][YY] <= 0));
 937
 938         if (epcPARRINELLORAHMAN == ir->epc && opts->bGenVel)
 939         {
 940             sprintf(warn_buf,
 941                     "You are generating velocities so I am assuming you "
 942                     "are equilibrating a system. You are using "
 943                     "%s pressure coupling, but this can be "
 944                     "unstable for equilibration. If your system crashes, try "
 945                     "equilibrating first with Berendsen pressure coupling. If "
 946                     "you are not equilibrating the system, you can probably "
 947                     "ignore this warning.",
 948                     epcoupl_names[ir->epc]);
 949             warning(wi,warn_buf);
 950         }
 951     }
 952
 953     if (EI_VV(ir->eI))
 954     {
 955         if (ir->epc > epcNO)
 956         {
 957             if ((ir->epc!=epcBERENDSEN) && (ir->epc!=epcMTTK))
 958             {
 959                 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.");
 960             }
 961         }
 962     }
 963
 964   /* ELECTROSTATICS */
 965   /* More checks are in triple check (grompp.c) */
 966
 967   if (ir->coulombtype == eelSWITCH) {
 968     sprintf(warn_buf,"coulombtype = %s is only for testing purposes and can lead to serious "
 969             "artifacts, advice: use coulombtype = %s",
 970             eel_names[ir->coulombtype],
 971             eel_names[eelRF_ZERO]);
 972     warning(wi,warn_buf);
 973   }
 974
 975   if (ir->epsilon_r!=1 && ir->implicit_solvent==eisGBSA) {
 976     sprintf(warn_buf,"epsilon-r = %g with GB implicit solvent, will use this value for inner dielectric",ir->epsilon_r);
 977     warning_note(wi,warn_buf);
 978   }
 979
 980   if (EEL_RF(ir->coulombtype) && ir->epsilon_rf==1 && ir->epsilon_r!=1) {
 981     sprintf(warn_buf,"epsilon-r = %g and epsilon-rf = 1 with reaction field, proceeding assuming old format and exchanging epsilon-r and epsilon-rf",ir->epsilon_r);
 982     warning(wi,warn_buf);
 983     ir->epsilon_rf = ir->epsilon_r;
 984     ir->epsilon_r  = 1.0;
 985   }
 986
 987   if (getenv("GALACTIC_DYNAMICS") == NULL) {
 988     sprintf(err_buf,"epsilon-r must be >= 0 instead of %g\n",ir->epsilon_r);
 989     CHECK(ir->epsilon_r < 0);
 990   }
 991
 992   if (EEL_RF(ir->coulombtype)) {
 993     /* reaction field (at the cut-off) */
 994
 995     if (ir->coulombtype == eelRF_ZERO) {
 996        sprintf(warn_buf,"With coulombtype = %s, epsilon-rf must be 0, assuming you meant epsilon_rf=0",
 997                eel_names[ir->coulombtype]);
 998         CHECK(ir->epsilon_rf != 0);
 999         ir->epsilon_rf = 0.0;
1000     }
1001
1002     sprintf(err_buf,"epsilon-rf must be >= epsilon-r");
1003     CHECK((ir->epsilon_rf < ir->epsilon_r && ir->epsilon_rf != 0) ||
1004           (ir->epsilon_r == 0));
1005     if (ir->epsilon_rf == ir->epsilon_r) {
1006       sprintf(warn_buf,"Using epsilon-rf = epsilon-r with %s does not make sense",
1007               eel_names[ir->coulombtype]);
1008       warning(wi,warn_buf);
1009     }
1010   }
1011   /* Allow rlist>rcoulomb for tabulated long range stuff. This just
1012    * means the interaction is zero outside rcoulomb, but it helps to
1013    * provide accurate energy conservation.
1014    */
1015   if (EEL_MIGHT_BE_ZERO_AT_CUTOFF(ir->coulombtype)) {
1016     if (EEL_SWITCHED(ir->coulombtype)) {
1017       sprintf(err_buf,
1018               "With coulombtype = %s rcoulomb_switch must be < rcoulomb. Or, better: Use the potential modifier options!",
1019               eel_names[ir->coulombtype]);
1020       CHECK(ir->rcoulomb_switch >= ir->rcoulomb);
1021     }
1022   } else if (ir->coulombtype == eelCUT || EEL_RF(ir->coulombtype)) {
1023       if (ir->cutoff_scheme == ecutsGROUP && ir->coulomb_modifier == eintmodNONE) {
1024           sprintf(err_buf,"With coulombtype = %s, rcoulomb should be >= rlist unless you use a potential modifier",
1025                   eel_names[ir->coulombtype]);
1026           CHECK(ir->rlist > ir->rcoulomb);
1027       }
1028   }
1029
1030   if(ir->coulombtype==eelSWITCH || ir->coulombtype==eelSHIFT ||
1031      ir->vdwtype==evdwSWITCH || ir->vdwtype==evdwSHIFT)
1032   {
1033       sprintf(warn_buf,
1034               "The switch/shift interaction settings are just for compatibility; you will get better"
1035               "performance from applying potential modifiers to your interactions!\n");
1036       warning_note(wi,warn_buf);
1037   }
1038
1039   if (EEL_FULL(ir->coulombtype))
1040   {
1041       if (ir->coulombtype==eelPMESWITCH || ir->coulombtype==eelPMEUSER ||
1042           ir->coulombtype==eelPMEUSERSWITCH)
1043       {
1044           sprintf(err_buf,"With coulombtype = %s, rcoulomb must be <= rlist",
1045                   eel_names[ir->coulombtype]);
1046           CHECK(ir->rcoulomb > ir->rlist);
1047       }
1048       else if (ir->cutoff_scheme == ecutsGROUP && ir->coulomb_modifier == eintmodNONE)
1049       {
1050           if (ir->coulombtype == eelPME || ir->coulombtype == eelP3M_AD)
1051           {
1052               sprintf(err_buf,
1053                       "With coulombtype = %s (without modifier), rcoulomb must be equal to rlist,\n"
1054                       "or rlistlong if nstcalclr=1. For optimal energy conservation,consider using\n"
1055                       "a potential modifier.",eel_names[ir->coulombtype]);
1056               if(ir->nstcalclr==1)
1057               {
1058                   CHECK(ir->rcoulomb != ir->rlist && ir->rcoulomb != ir->rlistlong);
1059               }
1060               else
1061               {
1062                   CHECK(ir->rcoulomb != ir->rlist);
1063               }
1064           }
1065       }
1066   }
1067
1068   if (EEL_PME(ir->coulombtype)) {
1069     if (ir->pme_order < 3) {
1070         warning_error(wi,"pme-order can not be smaller than 3");
1071     }
1072   }
1073
1074   if (ir->nwall==2 && EEL_FULL(ir->coulombtype)) {
1075     if (ir->ewald_geometry == eewg3D) {
1076       sprintf(warn_buf,"With pbc=%s you should use ewald-geometry=%s",
1077               epbc_names[ir->ePBC],eewg_names[eewg3DC]);
1078       warning(wi,warn_buf);
1079     }
1080     /* This check avoids extra pbc coding for exclusion corrections */
1081     sprintf(err_buf,"wall-ewald-zfac should be >= 2");
1082     CHECK(ir->wall_ewald_zfac < 2);
1083   }
1084
1085   if (EVDW_SWITCHED(ir->vdwtype)) {
1086     sprintf(err_buf,"With vdwtype = %s rvdw-switch must be < rvdw. Or, better - use a potential modifier.",
1087             evdw_names[ir->vdwtype]);
1088     CHECK(ir->rvdw_switch >= ir->rvdw);
1089   } else if (ir->vdwtype == evdwCUT) {
1090       if (ir->cutoff_scheme == ecutsGROUP && ir->vdw_modifier == eintmodNONE) {
1091           sprintf(err_buf,"With vdwtype = %s, rvdw must be >= rlist unless you use a potential modifier",evdw_names[ir->vdwtype]);
1092           CHECK(ir->rlist > ir->rvdw);
1093       }
1094   }
1095     if (ir->cutoff_scheme == ecutsGROUP)
1096     {
1097         if (EEL_IS_ZERO_AT_CUTOFF(ir->coulombtype)
1098             && (ir->rlistlong <= ir->rcoulomb))
1099         {
1100             sprintf(warn_buf,"For energy conservation with switch/shift potentials, %s should be 0.1 to 0.3 nm larger than rcoulomb.",
1101                     IR_TWINRANGE(*ir) ? "rlistlong" : "rlist");
1102             warning_note(wi,warn_buf);
1103         }
1104         if (EVDW_SWITCHED(ir->vdwtype) && (ir->rlistlong <= ir->rvdw))
1105         {
1106             sprintf(warn_buf,"For energy conservation with switch/shift potentials, %s should be 0.1 to 0.3 nm larger than rvdw.",
1107                     IR_TWINRANGE(*ir) ? "rlistlong" : "rlist");
1108             warning_note(wi,warn_buf);
1109         }
1110     }
1111
1112   if (ir->vdwtype == evdwUSER && ir->eDispCorr != edispcNO) {
1113       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.");
1114   }
1115
1116   if (ir->nstlist == -1) {
1117     sprintf(err_buf,"With nstlist=-1 rvdw and rcoulomb should be smaller than rlist to account for diffusion and possibly charge-group radii");
1118     CHECK(ir->rvdw >= ir->rlist || ir->rcoulomb >= ir->rlist);
1119   }
1120   sprintf(err_buf,"nstlist can not be smaller than -1");
1121   CHECK(ir->nstlist < -1);
1122
1123   if (ir->eI == eiLBFGS && (ir->coulombtype==eelCUT || ir->vdwtype==evdwCUT)
1124      && ir->rvdw != 0) {
1125     warning(wi,"For efficient BFGS minimization, use switch/shift/pme instead of cut-off.");
1126   }
1127
1128   if (ir->eI == eiLBFGS && ir->nbfgscorr <= 0) {
1129     warning(wi,"Using L-BFGS with nbfgscorr<=0 just gets you steepest descent.");
1130   }
1131
1132     /* ENERGY CONSERVATION */
1133     if (ir_NVE(ir) && ir->cutoff_scheme == ecutsGROUP)
1134     {
1135         if (!EVDW_MIGHT_BE_ZERO_AT_CUTOFF(ir->vdwtype) && ir->rvdw > 0 && ir->vdw_modifier == eintmodNONE)
1136         {
1137             sprintf(warn_buf,"You are using a cut-off for VdW interactions with NVE, for good energy conservation use vdwtype = %s (possibly with DispCorr)",
1138                     evdw_names[evdwSHIFT]);
1139             warning_note(wi,warn_buf);
1140         }
1141         if (!EEL_MIGHT_BE_ZERO_AT_CUTOFF(ir->coulombtype) && ir->rcoulomb > 0 && ir->coulomb_modifier == eintmodNONE)
1142         {
1143             sprintf(warn_buf,"You are using a cut-off for electrostatics with NVE, for good energy conservation use coulombtype = %s or %s",
1144                     eel_names[eelPMESWITCH],eel_names[eelRF_ZERO]);
1145             warning_note(wi,warn_buf);
1146         }
1147     }
1148
1149   /* IMPLICIT SOLVENT */
1150   if(ir->coulombtype==eelGB_NOTUSED)
1151   {
1152     ir->coulombtype=eelCUT;
1153     ir->implicit_solvent=eisGBSA;
1154     fprintf(stderr,"Note: Old option for generalized born electrostatics given:\n"
1155             "Changing coulombtype from \"generalized-born\" to \"cut-off\" and instead\n"
1156             "setting implicit-solvent value to \"GBSA\" in input section.\n");
1157   }
1158
1159   if(ir->sa_algorithm==esaSTILL)
1160   {
1161     sprintf(err_buf,"Still SA algorithm not available yet, use %s or %s instead\n",esa_names[esaAPPROX],esa_names[esaNO]);
1162     CHECK(ir->sa_algorithm == esaSTILL);
1163   }
1164
1165   if(ir->implicit_solvent==eisGBSA)
1166   {
1167     sprintf(err_buf,"With GBSA implicit solvent, rgbradii must be equal to rlist.");
1168     CHECK(ir->rgbradii != ir->rlist);
1169
1170     if(ir->coulombtype!=eelCUT)
1171           {
1172                   sprintf(err_buf,"With GBSA, coulombtype must be equal to %s\n",eel_names[eelCUT]);
1173                   CHECK(ir->coulombtype!=eelCUT);
1174           }
1175           if(ir->vdwtype!=evdwCUT)
1176           {
1177                   sprintf(err_buf,"With GBSA, vdw-type must be equal to %s\n",evdw_names[evdwCUT]);
1178                   CHECK(ir->vdwtype!=evdwCUT);
1179           }
1180     if(ir->nstgbradii<1)
1181     {
1182       sprintf(warn_buf,"Using GBSA with nstgbradii<1, setting nstgbradii=1");
1183       warning_note(wi,warn_buf);
1184       ir->nstgbradii=1;
1185     }
1186     if(ir->sa_algorithm==esaNO)
1187     {
1188       sprintf(warn_buf,"No SA (non-polar) calculation requested together with GB. Are you sure this is what you want?\n");
1189       warning_note(wi,warn_buf);
1190     }
1191     if(ir->sa_surface_tension<0 && ir->sa_algorithm!=esaNO)
1192     {
1193       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");
1194       warning_note(wi,warn_buf);
1195
1196       if(ir->gb_algorithm==egbSTILL)
1197       {
1198         ir->sa_surface_tension = 0.0049 * CAL2JOULE * 100;
1199       }
1200       else
1201       {
1202         ir->sa_surface_tension = 0.0054 * CAL2JOULE * 100;
1203       }
1204     }
1205     if(ir->sa_surface_tension==0 && ir->sa_algorithm!=esaNO)
1206     {
1207       sprintf(err_buf, "Surface tension set to 0 while SA-calculation requested\n");
1208       CHECK(ir->sa_surface_tension==0 && ir->sa_algorithm!=esaNO);
1209     }
1210
1211   }
1212
1213     if (ir->bAdress)
1214     {
1215         if (ir->cutoff_scheme != ecutsGROUP)
1216         {
1217             warning_error(wi,"AdresS simulation supports only cutoff-scheme=group");
1218         }
1219         if (!EI_SD(ir->eI))
1220         {
1221             warning_error(wi,"AdresS simulation supports only stochastic dynamics");
1222         }
1223         if (ir->epc != epcNO)
1224         {
1225             warning_error(wi,"AdresS simulation does not support pressure coupling");
1226         }
1227         if (EEL_FULL(ir->coulombtype))
1228         {
1229             warning_error(wi,"AdresS simulation does not support long-range electrostatics");
1230         }
1231     }
1232 }
1233
1234 /* count the number of text elemets separated by whitespace in a string.
1235     str = the input string
1236     maxptr = the maximum number of allowed elements
1237     ptr = the output array of pointers to the first character of each element
1238     returns: the number of elements. */
1239 int str_nelem(const char *str,int maxptr,char *ptr[])
1240 {
1241   int  np=0;
1242   char *copy0,*copy;
1243
1244   copy0=strdup(str);
1245   copy=copy0;
1246   ltrim(copy);
1247   while (*copy != '\0') {
1248     if (np >= maxptr)
1249       gmx_fatal(FARGS,"Too many groups on line: '%s' (max is %d)",
1250                   str,maxptr);
1251     if (ptr)
1252       ptr[np]=copy;
1253     np++;
1254     while ((*copy != '\0') && !isspace(*copy))
1255       copy++;
1256     if (*copy != '\0') {
1257       *copy='\0';
1258       copy++;
1259     }
1260     ltrim(copy);
1261   }
1262   if (ptr == NULL)
1263     sfree(copy0);
1264
1265   return np;
1266 }
1267
1268 /* interpret a number of doubles from a string and put them in an array,
1269    after allocating space for them.
1270    str = the input string
1271    n = the (pre-allocated) number of doubles read
1272    r = the output array of doubles. */
1273 static void parse_n_real(char *str,int *n,real **r)
1274 {
1275   char *ptr[MAXPTR];
1276   int  i;
1277
1278   *n = str_nelem(str,MAXPTR,ptr);
1279
1280   snew(*r,*n);
1281   for(i=0; i<*n; i++) {
1282     (*r)[i] = strtod(ptr[i],NULL);
1283   }
1284 }
1285
1286 static void do_fep_params(t_inputrec *ir, char fep_lambda[][STRLEN],char weights[STRLEN]) {
1287
1288     int i,j,max_n_lambda,nweights,nfep[efptNR];
1289     t_lambda *fep = ir->fepvals;
1290     t_expanded *expand = ir->expandedvals;
1291     real **count_fep_lambdas;
1292     gmx_bool bOneLambda = TRUE;
1293
1294     snew(count_fep_lambdas,efptNR);
1295
1296     /* FEP input processing */
1297     /* first, identify the number of lambda values for each type.
1298        All that are nonzero must have the same number */
1299
1300     for (i=0;i<efptNR;i++)
1301     {
1302         parse_n_real(fep_lambda[i],&(nfep[i]),&(count_fep_lambdas[i]));
1303     }
1304
1305     /* now, determine the number of components.  All must be either zero, or equal. */
1306
1307     max_n_lambda = 0;
1308     for (i=0;i<efptNR;i++)
1309     {
1310         if (nfep[i] > max_n_lambda) {
1311             max_n_lambda = nfep[i];  /* here's a nonzero one.  All of them
1312                                         must have the same number if its not zero.*/
1313             break;
1314         }
1315     }
1316
1317     for (i=0;i<efptNR;i++)
1318     {
1319         if (nfep[i] == 0)
1320         {
1321             ir->fepvals->separate_dvdl[i] = FALSE;
1322         }
1323         else if (nfep[i] == max_n_lambda)
1324         {
1325             if (i!=efptTEMPERATURE)  /* we treat this differently -- not really a reason to compute the derivative with
1326                                         respect to the temperature currently */
1327             {
1328                 ir->fepvals->separate_dvdl[i] = TRUE;
1329             }
1330         }
1331         else
1332         {
1333             gmx_fatal(FARGS,"Number of lambdas (%d) for FEP type %s not equal to number of other types (%d)",
1334                       nfep[i],efpt_names[i],max_n_lambda);
1335         }
1336     }
1337     /* we don't print out dhdl if the temperature is changing, since we can't correctly define dhdl in this case */
1338     ir->fepvals->separate_dvdl[efptTEMPERATURE] = FALSE;
1339
1340     /* the number of lambdas is the number we've read in, which is either zero
1341        or the same for all */
1342     fep->n_lambda = max_n_lambda;
1343
1344     /* allocate space for the array of lambda values */
1345     snew(fep->all_lambda,efptNR);
1346     /* if init_lambda is defined, we need to set lambda */
1347     if ((fep->init_lambda > 0) && (fep->n_lambda == 0))
1348     {
1349         ir->fepvals->separate_dvdl[efptFEP] = TRUE;
1350     }
1351     /* otherwise allocate the space for all of the lambdas, and transfer the data */
1352     for (i=0;i<efptNR;i++)
1353     {
1354         snew(fep->all_lambda[i],fep->n_lambda);
1355         if (nfep[i] > 0)  /* if it's zero, then the count_fep_lambda arrays
1356                              are zero */
1357         {
1358             for (j=0;j<fep->n_lambda;j++)
1359             {
1360                 fep->all_lambda[i][j] = (double)count_fep_lambdas[i][j];
1361             }
1362             sfree(count_fep_lambdas[i]);
1363         }
1364     }
1365     sfree(count_fep_lambdas);
1366
1367     /* "fep-vals" is either zero or the full number. If zero, we'll need to define fep-lambdas for internal
1368        bookkeeping -- for now, init_lambda */
1369
1370     if ((nfep[efptFEP] == 0) && (fep->init_lambda >= 0))
1371     {
1372         for (i=0;i<fep->n_lambda;i++)
1373         {
1374             fep->all_lambda[efptFEP][i] = fep->init_lambda;
1375         }
1376     }
1377
1378     /* check to see if only a single component lambda is defined, and soft core is defined.
1379        In this case, turn on coulomb soft core */
1380
1381     if (max_n_lambda == 0)
1382     {
1383         bOneLambda = TRUE;
1384     }
1385     else
1386     {
1387         for (i=0;i<efptNR;i++)
1388         {
1389             if ((nfep[i] != 0) && (i!=efptFEP))
1390             {
1391                 bOneLambda = FALSE;
1392             }
1393         }
1394     }
1395     if ((bOneLambda) && (fep->sc_alpha > 0))
1396     {
1397         fep->bScCoul = TRUE;
1398     }
1399
1400     /* Fill in the others with the efptFEP if they are not explicitly
1401        specified (i.e. nfep[i] == 0).  This means if fep is not defined,
1402        they are all zero. */
1403
1404     for (i=0;i<efptNR;i++)
1405     {
1406         if ((nfep[i] == 0) && (i!=efptFEP))
1407         {
1408             for (j=0;j<fep->n_lambda;j++)
1409             {
1410                 fep->all_lambda[i][j] = fep->all_lambda[efptFEP][j];
1411             }
1412         }
1413     }
1414
1415
1416     /* make it easier if sc_r_power = 48 by increasing it to the 4th power, to be in the right scale. */
1417     if (fep->sc_r_power == 48)
1418     {
1419         if (fep->sc_alpha > 0.1)
1420         {
1421             gmx_fatal(FARGS,"sc_alpha (%f) for sc_r_power = 48 should usually be between 0.001 and 0.004", fep->sc_alpha);
1422         }
1423     }
1424
1425     expand = ir->expandedvals;
1426     /* now read in the weights */
1427     parse_n_real(weights,&nweights,&(expand->init_lambda_weights));
1428     if (nweights == 0)
1429     {
1430         expand->bInit_weights = FALSE;
1431         snew(expand->init_lambda_weights,fep->n_lambda); /* initialize to zero */
1432     }
1433     else if (nweights != fep->n_lambda)
1434     {
1435         gmx_fatal(FARGS,"Number of weights (%d) is not equal to number of lambda values (%d)",
1436                   nweights,fep->n_lambda);
1437     }
1438     else
1439     {
1440         expand->bInit_weights = TRUE;
1441     }
1442     if ((expand->nstexpanded < 0) && (ir->efep != efepNO)) {
1443         expand->nstexpanded = fep->nstdhdl;
1444         /* if you don't specify nstexpanded when doing expanded ensemble free energy calcs, it is set to nstdhdl */
1445     }
1446     if ((expand->nstexpanded < 0) && ir->bSimTemp) {
1447         expand->nstexpanded = 2*(int)(ir->opts.tau_t[0]/ir->delta_t);
1448         /* if you don't specify nstexpanded when doing expanded ensemble simulated tempering, it is set to
1449            2*tau_t just to be careful so it's not to frequent  */
1450     }
1451 }
1452
1453
1454 static void do_simtemp_params(t_inputrec *ir) {
1455
1456     snew(ir->simtempvals->temperatures,ir->fepvals->n_lambda);
1457     GetSimTemps(ir->fepvals->n_lambda,ir->simtempvals,ir->fepvals->all_lambda[efptTEMPERATURE]);
1458
1459     return;
1460 }
1461
1462 static void do_wall_params(t_inputrec *ir,
1463                            char *wall_atomtype, char *wall_density,
1464                            t_gromppopts *opts)
1465 {
1466     int  nstr,i;
1467     char *names[MAXPTR];
1468     double dbl;
1469
1470     opts->wall_atomtype[0] = NULL;
1471     opts->wall_atomtype[1] = NULL;
1472
1473     ir->wall_atomtype[0] = -1;
1474     ir->wall_atomtype[1] = -1;
1475     ir->wall_density[0] = 0;
1476     ir->wall_density[1] = 0;
1477
1478     if (ir->nwall > 0)
1479     {
1480         nstr = str_nelem(wall_atomtype,MAXPTR,names);
1481         if (nstr != ir->nwall)
1482         {
1483             gmx_fatal(FARGS,"Expected %d elements for wall_atomtype, found %d",
1484                       ir->nwall,nstr);
1485         }
1486         for(i=0; i<ir->nwall; i++)
1487         {
1488             opts->wall_atomtype[i] = strdup(names[i]);
1489         }
1490
1491         if (ir->wall_type == ewt93 || ir->wall_type == ewt104) {
1492             nstr = str_nelem(wall_density,MAXPTR,names);
1493             if (nstr != ir->nwall)
1494             {
1495                 gmx_fatal(FARGS,"Expected %d elements for wall-density, found %d",ir->nwall,nstr);
1496             }
1497             for(i=0; i<ir->nwall; i++)
1498             {
1499                 sscanf(names[i],"%lf",&dbl);
1500                 if (dbl <= 0)
1501                 {
1502                     gmx_fatal(FARGS,"wall-density[%d] = %f\n",i,dbl);
1503                 }
1504                 ir->wall_density[i] = dbl;
1505             }
1506         }
1507     }
1508 }
1509
1510 static void add_wall_energrps(gmx_groups_t *groups,int nwall,t_symtab *symtab)
1511 {
1512   int  i;
1513   t_grps *grps;
1514   char str[STRLEN];
1515
1516   if (nwall > 0) {
1517     srenew(groups->grpname,groups->ngrpname+nwall);
1518     grps = &(groups->grps[egcENER]);
1519     srenew(grps->nm_ind,grps->nr+nwall);
1520     for(i=0; i<nwall; i++) {
1521       sprintf(str,"wall%d",i);
1522       groups->grpname[groups->ngrpname] = put_symtab(symtab,str);
1523       grps->nm_ind[grps->nr++] = groups->ngrpname++;
1524     }
1525   }
1526 }
1527
1528 void read_expandedparams(int *ninp_p,t_inpfile **inp_p,
1529                          t_expanded *expand,warninp_t wi)
1530 {
1531   int  ninp,nerror=0;
1532   t_inpfile *inp;
1533
1534   ninp   = *ninp_p;
1535   inp    = *inp_p;
1536
1537   /* read expanded ensemble parameters */
1538   CCTYPE ("expanded ensemble variables");
1539   ITYPE ("nstexpanded",expand->nstexpanded,-1);
1540   EETYPE("lmc-stats", expand->elamstats, elamstats_names);
1541   EETYPE("lmc-move", expand->elmcmove, elmcmove_names);
1542   EETYPE("lmc-weights-equil",expand->elmceq,elmceq_names);
1543   ITYPE ("weight-equil-number-all-lambda",expand->equil_n_at_lam,-1);
1544   ITYPE ("weight-equil-number-samples",expand->equil_samples,-1);
1545   ITYPE ("weight-equil-number-steps",expand->equil_steps,-1);
1546   RTYPE ("weight-equil-wl-delta",expand->equil_wl_delta,-1);
1547   RTYPE ("weight-equil-count-ratio",expand->equil_ratio,-1);
1548   CCTYPE("Seed for Monte Carlo in lambda space");
1549   ITYPE ("lmc-seed",expand->lmc_seed,-1);
1550   RTYPE ("mc-temperature",expand->mc_temp,-1);
1551   ITYPE ("lmc-repeats",expand->lmc_repeats,1);
1552   ITYPE ("lmc-gibbsdelta",expand->gibbsdeltalam,-1);
1553   ITYPE ("lmc-forced-nstart",expand->lmc_forced_nstart,0);
1554   EETYPE("symmetrized-transition-matrix", expand->bSymmetrizedTMatrix, yesno_names);
1555   ITYPE("nst-transition-matrix", expand->nstTij, -1);
1556   ITYPE ("mininum-var-min",expand->minvarmin, 100); /*default is reasonable */
1557   ITYPE ("weight-c-range",expand->c_range, 0); /* default is just C=0 */
1558   RTYPE ("wl-scale",expand->wl_scale,0.8);
1559   RTYPE ("wl-ratio",expand->wl_ratio,0.8);
1560   RTYPE ("init-wl-delta",expand->init_wl_delta,1.0);
1561   EETYPE("wl-oneovert",expand->bWLoneovert,yesno_names);
1562
1563   *ninp_p   = ninp;
1564   *inp_p    = inp;
1565
1566   return;
1567 }
1568
1569 void get_ir(const char *mdparin,const char *mdparout,
1570             t_inputrec *ir,t_gromppopts *opts,
1571             warninp_t wi)
1572 {
1573   char      *dumstr[2];
1574   double    dumdub[2][6];
1575   t_inpfile *inp;
1576   const char *tmp;
1577   int       i,j,m,ninp;
1578   char      warn_buf[STRLEN];
1579   t_lambda  *fep = ir->fepvals;
1580   t_expanded *expand = ir->expandedvals;
1581
1582   inp = read_inpfile(mdparin, &ninp, NULL, wi);
1583
1584   snew(dumstr[0],STRLEN);
1585   snew(dumstr[1],STRLEN);
1586
1587   /* remove the following deprecated commands */
1588   REM_TYPE("title");
1589   REM_TYPE("cpp");
1590   REM_TYPE("domain-decomposition");
1591   REM_TYPE("andersen-seed");
1592   REM_TYPE("dihre");
1593   REM_TYPE("dihre-fc");
1594   REM_TYPE("dihre-tau");
1595   REM_TYPE("nstdihreout");
1596   REM_TYPE("nstcheckpoint");
1597
1598   /* replace the following commands with the clearer new versions*/
1599   REPL_TYPE("unconstrained-start","continuation");
1600   REPL_TYPE("foreign-lambda","fep-lambdas");
1601
1602   CCTYPE ("VARIOUS PREPROCESSING OPTIONS");
1603   CTYPE ("Preprocessor information: use cpp syntax.");
1604   CTYPE ("e.g.: -I/home/joe/doe -I/home/mary/roe");
1605   STYPE ("include",     opts->include,  NULL);
1606   CTYPE ("e.g.: -DPOSRES -DFLEXIBLE (note these variable names are case sensitive)");
1607   STYPE ("define",      opts->define,   NULL);
1608
1609   CCTYPE ("RUN CONTROL PARAMETERS");
1610   EETYPE("integrator",  ir->eI,         ei_names);
1611   CTYPE ("Start time and timestep in ps");
1612   RTYPE ("tinit",       ir->init_t,     0.0);
1613   RTYPE ("dt",          ir->delta_t,    0.001);
1614   STEPTYPE ("nsteps",   ir->nsteps,     0);
1615   CTYPE ("For exact run continuation or redoing part of a run");
1616   STEPTYPE ("init-step",ir->init_step,  0);
1617   CTYPE ("Part index is updated automatically on checkpointing (keeps files separate)");
1618   ITYPE ("simulation-part", ir->simulation_part, 1);
1619   CTYPE ("mode for center of mass motion removal");
1620   EETYPE("comm-mode",   ir->comm_mode,  ecm_names);
1621   CTYPE ("number of steps for center of mass motion removal");
1622   ITYPE ("nstcomm",     ir->nstcomm,    100);
1623   CTYPE ("group(s) for center of mass motion removal");
1624   STYPE ("comm-grps",   vcm,            NULL);
1625
1626   CCTYPE ("LANGEVIN DYNAMICS OPTIONS");
1627   CTYPE ("Friction coefficient (amu/ps) and random seed");
1628   RTYPE ("bd-fric",     ir->bd_fric,    0.0);
1629   ITYPE ("ld-seed",     ir->ld_seed,    1993);
1630
1631   /* Em stuff */
1632   CCTYPE ("ENERGY MINIMIZATION OPTIONS");
1633   CTYPE ("Force tolerance and initial step-size");
1634   RTYPE ("emtol",       ir->em_tol,     10.0);
1635   RTYPE ("emstep",      ir->em_stepsize,0.01);
1636   CTYPE ("Max number of iterations in relax-shells");
1637   ITYPE ("niter",       ir->niter,      20);
1638   CTYPE ("Step size (ps^2) for minimization of flexible constraints");
1639   RTYPE ("fcstep",      ir->fc_stepsize, 0);
1640   CTYPE ("Frequency of steepest descents steps when doing CG");
1641   ITYPE ("nstcgsteep",  ir->nstcgsteep, 1000);
1642   ITYPE ("nbfgscorr",   ir->nbfgscorr,  10);
1643
1644   CCTYPE ("TEST PARTICLE INSERTION OPTIONS");
1645   RTYPE ("rtpi",        ir->rtpi,       0.05);
1646
1647   /* Output options */
1648   CCTYPE ("OUTPUT CONTROL OPTIONS");
1649   CTYPE ("Output frequency for coords (x), velocities (v) and forces (f)");
1650   ITYPE ("nstxout",     ir->nstxout,    0);
1651   ITYPE ("nstvout",     ir->nstvout,    0);
1652   ITYPE ("nstfout",     ir->nstfout,    0);
1653   ir->nstcheckpoint = 1000;
1654   CTYPE ("Output frequency for energies to log file and energy file");
1655   ITYPE ("nstlog",      ir->nstlog,     1000);
1656   ITYPE ("nstcalcenergy",ir->nstcalcenergy,     100);
1657   ITYPE ("nstenergy",   ir->nstenergy,  1000);
1658   CTYPE ("Output frequency and precision for .xtc file");
1659   ITYPE ("nstxtcout",   ir->nstxtcout,  0);
1660   RTYPE ("xtc-precision",ir->xtcprec,   1000.0);
1661   CTYPE ("This selects the subset of atoms for the .xtc file. You can");
1662   CTYPE ("select multiple groups. By default all atoms will be written.");
1663   STYPE ("xtc-grps",    xtc_grps,       NULL);
1664   CTYPE ("Selection of energy groups");
1665   STYPE ("energygrps",  energy,         NULL);
1666
1667   /* Neighbor searching */
1668   CCTYPE ("NEIGHBORSEARCHING PARAMETERS");
1669   CTYPE ("cut-off scheme (group: using charge groups, Verlet: particle based cut-offs)");
1670   EETYPE("cutoff-scheme",     ir->cutoff_scheme,    ecutscheme_names);
1671   CTYPE ("nblist update frequency");
1672   ITYPE ("nstlist",     ir->nstlist,    10);
1673   CTYPE ("ns algorithm (simple or grid)");
1674   EETYPE("ns-type",     ir->ns_type,    ens_names);
1675   /* set ndelta to the optimal value of 2 */
1676   ir->ndelta = 2;
1677   CTYPE ("Periodic boundary conditions: xyz, no, xy");
1678   EETYPE("pbc",         ir->ePBC,       epbc_names);
1679   EETYPE("periodic-molecules", ir->bPeriodicMols, yesno_names);
1680   CTYPE ("Allowed energy drift due to the Verlet buffer in kJ/mol/ps per atom,");
1681   CTYPE ("a value of -1 means: use rlist");
1682   RTYPE("verlet-buffer-drift", ir->verletbuf_drift,    0.005);
1683   CTYPE ("nblist cut-off");
1684   RTYPE ("rlist",       ir->rlist,      1.0);
1685   CTYPE ("long-range cut-off for switched potentials");
1686   RTYPE ("rlistlong",   ir->rlistlong,  -1);
1687   ITYPE ("nstcalclr",   ir->nstcalclr,  -1);
1688
1689   /* Electrostatics */
1690   CCTYPE ("OPTIONS FOR ELECTROSTATICS AND VDW");
1691   CTYPE ("Method for doing electrostatics");
1692   EETYPE("coulombtype", ir->coulombtype,    eel_names);
1693   EETYPE("coulomb-modifier",    ir->coulomb_modifier,    eintmod_names);
1694   CTYPE ("cut-off lengths");
1695   RTYPE ("rcoulomb-switch",     ir->rcoulomb_switch,    0.0);
1696   RTYPE ("rcoulomb",    ir->rcoulomb,   1.0);
1697   CTYPE ("Relative dielectric constant for the medium and the reaction field");
1698   RTYPE ("epsilon-r",   ir->epsilon_r,  1.0);
1699   RTYPE ("epsilon-rf",  ir->epsilon_rf, 0.0);
1700   CTYPE ("Method for doing Van der Waals");
1701   EETYPE("vdw-type",    ir->vdwtype,    evdw_names);
1702   EETYPE("vdw-modifier",        ir->vdw_modifier,    eintmod_names);
1703   CTYPE ("cut-off lengths");
1704   RTYPE ("rvdw-switch", ir->rvdw_switch,        0.0);
1705   RTYPE ("rvdw",        ir->rvdw,       1.0);
1706   CTYPE ("Apply long range dispersion corrections for Energy and Pressure");
1707   EETYPE("DispCorr",    ir->eDispCorr,  edispc_names);
1708   CTYPE ("Extension of the potential lookup tables beyond the cut-off");
1709   RTYPE ("table-extension", ir->tabext, 1.0);
1710   CTYPE ("Separate tables between energy group pairs");
1711   STYPE ("energygrp-table", egptable,   NULL);
1712   CTYPE ("Spacing for the PME/PPPM FFT grid");
1713   RTYPE ("fourierspacing", ir->fourier_spacing,0.12);
1714   CTYPE ("FFT grid size, when a value is 0 fourierspacing will be used");
1715   ITYPE ("fourier-nx",  ir->nkx,         0);
1716   ITYPE ("fourier-ny",  ir->nky,         0);
1717   ITYPE ("fourier-nz",  ir->nkz,         0);
1718   CTYPE ("EWALD/PME/PPPM parameters");
1719   ITYPE ("pme-order",   ir->pme_order,   4);
1720   RTYPE ("ewald-rtol",  ir->ewald_rtol, 0.00001);
1721   EETYPE("ewald-geometry", ir->ewald_geometry, eewg_names);
1722   RTYPE ("epsilon-surface", ir->epsilon_surface, 0.0);
1723   EETYPE("optimize-fft",ir->bOptFFT,  yesno_names);
1724
1725   CCTYPE("IMPLICIT SOLVENT ALGORITHM");
1726   EETYPE("implicit-solvent", ir->implicit_solvent, eis_names);
1727
1728   CCTYPE ("GENERALIZED BORN ELECTROSTATICS");
1729   CTYPE ("Algorithm for calculating Born radii");
1730   EETYPE("gb-algorithm", ir->gb_algorithm, egb_names);
1731   CTYPE ("Frequency of calculating the Born radii inside rlist");
1732   ITYPE ("nstgbradii", ir->nstgbradii, 1);
1733   CTYPE ("Cutoff for Born radii calculation; the contribution from atoms");
1734   CTYPE ("between rlist and rgbradii is updated every nstlist steps");
1735   RTYPE ("rgbradii",  ir->rgbradii, 1.0);
1736   CTYPE ("Dielectric coefficient of the implicit solvent");
1737   RTYPE ("gb-epsilon-solvent",ir->gb_epsilon_solvent, 80.0);
1738   CTYPE ("Salt concentration in M for Generalized Born models");
1739   RTYPE ("gb-saltconc",  ir->gb_saltconc, 0.0);
1740   CTYPE ("Scaling factors used in the OBC GB model. Default values are OBC(II)");
1741   RTYPE ("gb-obc-alpha", ir->gb_obc_alpha, 1.0);
1742   RTYPE ("gb-obc-beta", ir->gb_obc_beta, 0.8);
1743   RTYPE ("gb-obc-gamma", ir->gb_obc_gamma, 4.85);
1744   RTYPE ("gb-dielectric-offset", ir->gb_dielectric_offset, 0.009);
1745   EETYPE("sa-algorithm", ir->sa_algorithm, esa_names);
1746   CTYPE ("Surface tension (kJ/mol/nm^2) for the SA (nonpolar surface) part of GBSA");
1747   CTYPE ("The value -1 will set default value for Still/HCT/OBC GB-models.");
1748   RTYPE ("sa-surface-tension", ir->sa_surface_tension, -1);
1749
1750   /* Coupling stuff */
1751   CCTYPE ("OPTIONS FOR WEAK COUPLING ALGORITHMS");
1752   CTYPE ("Temperature coupling");
1753   EETYPE("tcoupl",      ir->etc,        etcoupl_names);
1754   ITYPE ("nsttcouple", ir->nsttcouple,  -1);
1755   ITYPE("nh-chain-length",     ir->opts.nhchainlength, NHCHAINLENGTH);
1756   EETYPE("print-nose-hoover-chain-variables", ir->bPrintNHChains, yesno_names);
1757   CTYPE ("Groups to couple separately");
1758   STYPE ("tc-grps",     tcgrps,         NULL);
1759   CTYPE ("Time constant (ps) and reference temperature (K)");
1760   STYPE ("tau-t",       tau_t,          NULL);
1761   STYPE ("ref-t",       ref_t,          NULL);
1762   CTYPE ("pressure coupling");
1763   EETYPE("pcoupl",      ir->epc,        epcoupl_names);
1764   EETYPE("pcoupltype",  ir->epct,       epcoupltype_names);
1765   ITYPE ("nstpcouple", ir->nstpcouple,  -1);
1766   CTYPE ("Time constant (ps), compressibility (1/bar) and reference P (bar)");
1767   RTYPE ("tau-p",       ir->tau_p,      1.0);
1768   STYPE ("compressibility",     dumstr[0],      NULL);
1769   STYPE ("ref-p",       dumstr[1],      NULL);
1770   CTYPE ("Scaling of reference coordinates, No, All or COM");
1771   EETYPE ("refcoord-scaling",ir->refcoord_scaling,erefscaling_names);
1772
1773   /* QMMM */
1774   CCTYPE ("OPTIONS FOR QMMM calculations");
1775   EETYPE("QMMM", ir->bQMMM, yesno_names);
1776   CTYPE ("Groups treated Quantum Mechanically");
1777   STYPE ("QMMM-grps",  QMMM,          NULL);
1778   CTYPE ("QM method");
1779   STYPE("QMmethod",     QMmethod, NULL);
1780   CTYPE ("QMMM scheme");
1781   EETYPE("QMMMscheme",  ir->QMMMscheme,    eQMMMscheme_names);
1782   CTYPE ("QM basisset");
1783   STYPE("QMbasis",      QMbasis, NULL);
1784   CTYPE ("QM charge");
1785   STYPE ("QMcharge",    QMcharge,NULL);
1786   CTYPE ("QM multiplicity");
1787   STYPE ("QMmult",      QMmult,NULL);
1788   CTYPE ("Surface Hopping");
1789   STYPE ("SH",          bSH, NULL);
1790   CTYPE ("CAS space options");
1791   STYPE ("CASorbitals",      CASorbitals,   NULL);
1792   STYPE ("CASelectrons",     CASelectrons,  NULL);
1793   STYPE ("SAon", SAon, NULL);
1794   STYPE ("SAoff",SAoff,NULL);
1795   STYPE ("SAsteps",  SAsteps, NULL);
1796   CTYPE ("Scale factor for MM charges");
1797   RTYPE ("MMChargeScaleFactor", ir->scalefactor, 1.0);
1798   CTYPE ("Optimization of QM subsystem");
1799   STYPE ("bOPT",          bOPT, NULL);
1800   STYPE ("bTS",          bTS, NULL);
1801
1802   /* Simulated annealing */
1803   CCTYPE("SIMULATED ANNEALING");
1804   CTYPE ("Type of annealing for each temperature group (no/single/periodic)");
1805   STYPE ("annealing",   anneal,      NULL);
1806   CTYPE ("Number of time points to use for specifying annealing in each group");
1807   STYPE ("annealing-npoints", anneal_npoints, NULL);
1808   CTYPE ("List of times at the annealing points for each group");
1809   STYPE ("annealing-time",       anneal_time,       NULL);
1810   CTYPE ("Temp. at each annealing point, for each group.");
1811   STYPE ("annealing-temp",  anneal_temp,  NULL);
1812
1813   /* Startup run */
1814   CCTYPE ("GENERATE VELOCITIES FOR STARTUP RUN");
1815   EETYPE("gen-vel",     opts->bGenVel,  yesno_names);
1816   RTYPE ("gen-temp",    opts->tempi,    300.0);
1817   ITYPE ("gen-seed",    opts->seed,     173529);
1818
1819   /* Shake stuff */
1820   CCTYPE ("OPTIONS FOR BONDS");
1821   EETYPE("constraints", opts->nshake,   constraints);
1822   CTYPE ("Type of constraint algorithm");
1823   EETYPE("constraint-algorithm",  ir->eConstrAlg, econstr_names);
1824   CTYPE ("Do not constrain the start configuration");
1825   EETYPE("continuation", ir->bContinuation, yesno_names);
1826   CTYPE ("Use successive overrelaxation to reduce the number of shake iterations");
1827   EETYPE("Shake-SOR", ir->bShakeSOR, yesno_names);
1828   CTYPE ("Relative tolerance of shake");
1829   RTYPE ("shake-tol", ir->shake_tol, 0.0001);
1830   CTYPE ("Highest order in the expansion of the constraint coupling matrix");
1831   ITYPE ("lincs-order", ir->nProjOrder, 4);
1832   CTYPE ("Number of iterations in the final step of LINCS. 1 is fine for");
1833   CTYPE ("normal simulations, but use 2 to conserve energy in NVE runs.");
1834   CTYPE ("For energy minimization with constraints it should be 4 to 8.");
1835   ITYPE ("lincs-iter", ir->nLincsIter, 1);
1836   CTYPE ("Lincs will write a warning to the stderr if in one step a bond");
1837   CTYPE ("rotates over more degrees than");
1838   RTYPE ("lincs-warnangle", ir->LincsWarnAngle, 30.0);
1839   CTYPE ("Convert harmonic bonds to morse potentials");
1840   EETYPE("morse",       opts->bMorse,yesno_names);
1841
1842   /* Energy group exclusions */
1843   CCTYPE ("ENERGY GROUP EXCLUSIONS");
1844   CTYPE ("Pairs of energy groups for which all non-bonded interactions are excluded");
1845   STYPE ("energygrp-excl", egpexcl,     NULL);
1846
1847   /* Walls */
1848   CCTYPE ("WALLS");
1849   CTYPE ("Number of walls, type, atom types, densities and box-z scale factor for Ewald");
1850   ITYPE ("nwall", ir->nwall, 0);
1851   EETYPE("wall-type",     ir->wall_type,   ewt_names);
1852   RTYPE ("wall-r-linpot", ir->wall_r_linpot, -1);
1853   STYPE ("wall-atomtype", wall_atomtype, NULL);
1854   STYPE ("wall-density",  wall_density,  NULL);
1855   RTYPE ("wall-ewald-zfac", ir->wall_ewald_zfac, 3);
1856
1857   /* COM pulling */
1858   CCTYPE("COM PULLING");
1859   CTYPE("Pull type: no, umbrella, constraint or constant-force");
1860   EETYPE("pull",          ir->ePull, epull_names);
1861   if (ir->ePull != epullNO) {
1862     snew(ir->pull,1);
1863     pull_grp = read_pullparams(&ninp,&inp,ir->pull,&opts->pull_start,wi);
1864   }
1865
1866   /* Enforced rotation */
1867   CCTYPE("ENFORCED ROTATION");
1868   CTYPE("Enforced rotation: No or Yes");
1869   EETYPE("rotation",       ir->bRot, yesno_names);
1870   if (ir->bRot) {
1871     snew(ir->rot,1);
1872     rot_grp = read_rotparams(&ninp,&inp,ir->rot,wi);
1873   }
1874
1875   /* Refinement */
1876   CCTYPE("NMR refinement stuff");
1877   CTYPE ("Distance restraints type: No, Simple or Ensemble");
1878   EETYPE("disre",       ir->eDisre,     edisre_names);
1879   CTYPE ("Force weighting of pairs in one distance restraint: Conservative or Equal");
1880   EETYPE("disre-weighting", ir->eDisreWeighting, edisreweighting_names);
1881   CTYPE ("Use sqrt of the time averaged times the instantaneous violation");
1882   EETYPE("disre-mixed", ir->bDisreMixed, yesno_names);
1883   RTYPE ("disre-fc",    ir->dr_fc,      1000.0);
1884   RTYPE ("disre-tau",   ir->dr_tau,     0.0);
1885   CTYPE ("Output frequency for pair distances to energy file");
1886   ITYPE ("nstdisreout", ir->nstdisreout, 100);
1887   CTYPE ("Orientation restraints: No or Yes");
1888   EETYPE("orire",       opts->bOrire,   yesno_names);
1889   CTYPE ("Orientation restraints force constant and tau for time averaging");
1890   RTYPE ("orire-fc",    ir->orires_fc,  0.0);
1891   RTYPE ("orire-tau",   ir->orires_tau, 0.0);
1892   STYPE ("orire-fitgrp",orirefitgrp,    NULL);
1893   CTYPE ("Output frequency for trace(SD) and S to energy file");
1894   ITYPE ("nstorireout", ir->nstorireout, 100);
1895
1896   /* free energy variables */
1897   CCTYPE ("Free energy variables");
1898   EETYPE("free-energy", ir->efep, efep_names);
1899   STYPE ("couple-moltype",  couple_moltype,  NULL);
1900   EETYPE("couple-lambda0", opts->couple_lam0, couple_lam);
1901   EETYPE("couple-lambda1", opts->couple_lam1, couple_lam);
1902   EETYPE("couple-intramol", opts->bCoupleIntra, yesno_names);
1903
1904   RTYPE ("init-lambda", fep->init_lambda,-1); /* start with -1 so
1905                                                  we can recognize if
1906                                                  it was not entered */
1907   ITYPE ("init-lambda-state", fep->init_fep_state,-1);
1908   RTYPE ("delta-lambda",fep->delta_lambda,0.0);
1909   ITYPE ("nstdhdl",fep->nstdhdl, 50);
1910   STYPE ("fep-lambdas", fep_lambda[efptFEP], NULL);
1911   STYPE ("mass-lambdas", fep_lambda[efptMASS], NULL);
1912   STYPE ("coul-lambdas", fep_lambda[efptCOUL], NULL);
1913   STYPE ("vdw-lambdas", fep_lambda[efptVDW], NULL);
1914   STYPE ("bonded-lambdas", fep_lambda[efptBONDED], NULL);
1915   STYPE ("restraint-lambdas", fep_lambda[efptRESTRAINT], NULL);
1916   STYPE ("temperature-lambdas", fep_lambda[efptTEMPERATURE], NULL);
1917   ITYPE ("calc-lambda-neighbors",fep->lambda_neighbors, 1);
1918   STYPE ("init-lambda-weights",lambda_weights,NULL);
1919   EETYPE("dhdl-print-energy", fep->bPrintEnergy, yesno_names);
1920   RTYPE ("sc-alpha",fep->sc_alpha,0.0);
1921   ITYPE ("sc-power",fep->sc_power,1);
1922   RTYPE ("sc-r-power",fep->sc_r_power,6.0);
1923   RTYPE ("sc-sigma",fep->sc_sigma,0.3);
1924   EETYPE("sc-coul",fep->bScCoul,yesno_names);
1925   ITYPE ("dh_hist_size", fep->dh_hist_size, 0);
1926   RTYPE ("dh_hist_spacing", fep->dh_hist_spacing, 0.1);
1927   EETYPE("separate-dhdl-file", fep->separate_dhdl_file,
1928                                separate_dhdl_file_names);
1929   EETYPE("dhdl-derivatives", fep->dhdl_derivatives, dhdl_derivatives_names);
1930   ITYPE ("dh_hist_size", fep->dh_hist_size, 0);
1931   RTYPE ("dh_hist_spacing", fep->dh_hist_spacing, 0.1);
1932
1933   /* Non-equilibrium MD stuff */
1934   CCTYPE("Non-equilibrium MD stuff");
1935   STYPE ("acc-grps",    accgrps,        NULL);
1936   STYPE ("accelerate",  acc,            NULL);
1937   STYPE ("freezegrps",  freeze,         NULL);
1938   STYPE ("freezedim",   frdim,          NULL);
1939   RTYPE ("cos-acceleration", ir->cos_accel, 0);
1940   STYPE ("deform",      deform,         NULL);
1941
1942   /* simulated tempering variables */
1943   CCTYPE("simulated tempering variables");
1944   EETYPE("simulated-tempering",ir->bSimTemp,yesno_names);
1945   EETYPE("simulated-tempering-scaling",ir->simtempvals->eSimTempScale,esimtemp_names);
1946   RTYPE("sim-temp-low",ir->simtempvals->simtemp_low,300.0);
1947   RTYPE("sim-temp-high",ir->simtempvals->simtemp_high,300.0);
1948
1949   /* expanded ensemble variables */
1950   if (ir->efep==efepEXPANDED || ir->bSimTemp)
1951   {
1952       read_expandedparams(&ninp,&inp,expand,wi);
1953   }
1954
1955   /* Electric fields */
1956   CCTYPE("Electric fields");
1957   CTYPE ("Format is number of terms (int) and for all terms an amplitude (real)");
1958   CTYPE ("and a phase angle (real)");
1959   STYPE ("E-x",         efield_x,       NULL);
1960   STYPE ("E-xt",        efield_xt,      NULL);
1961   STYPE ("E-y",         efield_y,       NULL);
1962   STYPE ("E-yt",        efield_yt,      NULL);
1963   STYPE ("E-z",         efield_z,       NULL);
1964   STYPE ("E-zt",        efield_zt,      NULL);
1965
1966   /* AdResS defined thingies */
1967   CCTYPE ("AdResS parameters");
1968   EETYPE("adress",       ir->bAdress, yesno_names);
1969   if (ir->bAdress) {
1970     snew(ir->adress,1);
1971     read_adressparams(&ninp,&inp,ir->adress,wi);
1972   }
1973
1974   /* User defined thingies */
1975   CCTYPE ("User defined thingies");
1976   STYPE ("user1-grps",  user1,          NULL);
1977   STYPE ("user2-grps",  user2,          NULL);
1978   ITYPE ("userint1",    ir->userint1,   0);
1979   ITYPE ("userint2",    ir->userint2,   0);
1980   ITYPE ("userint3",    ir->userint3,   0);
1981   ITYPE ("userint4",    ir->userint4,   0);
1982   RTYPE ("userreal1",   ir->userreal1,  0);
1983   RTYPE ("userreal2",   ir->userreal2,  0);
1984   RTYPE ("userreal3",   ir->userreal3,  0);
1985   RTYPE ("userreal4",   ir->userreal4,  0);
1986 #undef CTYPE
1987
1988   write_inpfile(mdparout,ninp,inp,FALSE,wi);
1989   for (i=0; (i<ninp); i++) {
1990     sfree(inp[i].name);
1991     sfree(inp[i].value);
1992   }
1993   sfree(inp);
1994
1995   /* Process options if necessary */
1996   for(m=0; m<2; m++) {
1997     for(i=0; i<2*DIM; i++)
1998       dumdub[m][i]=0.0;
1999     if(ir->epc) {
2000       switch (ir->epct) {
2001       case epctISOTROPIC:
2002         if (sscanf(dumstr[m],"%lf",&(dumdub[m][XX]))!=1) {
2003         warning_error(wi,"Pressure coupling not enough values (I need 1)");
2004         }
2005         dumdub[m][YY]=dumdub[m][ZZ]=dumdub[m][XX];
2006         break;
2007       case epctSEMIISOTROPIC:
2008       case epctSURFACETENSION:
2009         if (sscanf(dumstr[m],"%lf%lf",
2010                    &(dumdub[m][XX]),&(dumdub[m][ZZ]))!=2) {
2011         warning_error(wi,"Pressure coupling not enough values (I need 2)");
2012         }
2013         dumdub[m][YY]=dumdub[m][XX];
2014         break;
2015       case epctANISOTROPIC:
2016         if (sscanf(dumstr[m],"%lf%lf%lf%lf%lf%lf",
2017                    &(dumdub[m][XX]),&(dumdub[m][YY]),&(dumdub[m][ZZ]),
2018                    &(dumdub[m][3]),&(dumdub[m][4]),&(dumdub[m][5]))!=6) {
2019         warning_error(wi,"Pressure coupling not enough values (I need 6)");
2020         }
2021         break;
2022       default:
2023         gmx_fatal(FARGS,"Pressure coupling type %s not implemented yet",
2024                     epcoupltype_names[ir->epct]);
2025       }
2026     }
2027   }
2028   clear_mat(ir->ref_p);
2029   clear_mat(ir->compress);
2030   for(i=0; i<DIM; i++) {
2031     ir->ref_p[i][i]    = dumdub[1][i];
2032     ir->compress[i][i] = dumdub[0][i];
2033   }
2034   if (ir->epct == epctANISOTROPIC) {
2035     ir->ref_p[XX][YY] = dumdub[1][3];
2036     ir->ref_p[XX][ZZ] = dumdub[1][4];
2037     ir->ref_p[YY][ZZ] = dumdub[1][5];
2038     if (ir->ref_p[XX][YY]!=0 && ir->ref_p[XX][ZZ]!=0 && ir->ref_p[YY][ZZ]!=0) {
2039       warning(wi,"All off-diagonal reference pressures are non-zero. Are you sure you want to apply a threefold shear stress?\n");
2040     }
2041     ir->compress[XX][YY] = dumdub[0][3];
2042     ir->compress[XX][ZZ] = dumdub[0][4];
2043     ir->compress[YY][ZZ] = dumdub[0][5];
2044     for(i=0; i<DIM; i++) {
2045       for(m=0; m<i; m++) {
2046         ir->ref_p[i][m] = ir->ref_p[m][i];
2047         ir->compress[i][m] = ir->compress[m][i];
2048       }
2049     }
2050   }
2051
2052   if (ir->comm_mode == ecmNO)
2053     ir->nstcomm = 0;
2054
2055   opts->couple_moltype = NULL;
2056   if (strlen(couple_moltype) > 0)
2057   {
2058       if (ir->efep != efepNO)
2059       {
2060           opts->couple_moltype = strdup(couple_moltype);
2061           if (opts->couple_lam0 == opts->couple_lam1)
2062           {
2063               warning(wi,"The lambda=0 and lambda=1 states for coupling are identical");
2064           }
2065           if (ir->eI == eiMD && (opts->couple_lam0 == ecouplamNONE ||
2066                                  opts->couple_lam1 == ecouplamNONE))
2067           {
2068               warning(wi,"For proper sampling of the (nearly) decoupled state, stochastic dynamics should be used");
2069           }
2070       }
2071       else
2072       {
2073           warning(wi,"Can not couple a molecule with free_energy = no");
2074       }
2075   }
2076   /* FREE ENERGY AND EXPANDED ENSEMBLE OPTIONS */
2077   if (ir->efep != efepNO) {
2078       if (fep->delta_lambda > 0) {
2079           ir->efep = efepSLOWGROWTH;
2080       }
2081   }
2082
2083   if (ir->bSimTemp) {
2084       fep->bPrintEnergy = TRUE;
2085       /* always print out the energy to dhdl if we are doing expanded ensemble, since we need the total energy
2086          if the temperature is changing. */
2087   }
2088
2089   if ((ir->efep != efepNO) || ir->bSimTemp)
2090   {
2091       ir->bExpanded = FALSE;
2092       if ((ir->efep == efepEXPANDED) || ir->bSimTemp)
2093       {
2094           ir->bExpanded = TRUE;
2095       }
2096       do_fep_params(ir,fep_lambda,lambda_weights);
2097       if (ir->bSimTemp) { /* done after fep params */
2098           do_simtemp_params(ir);
2099       }
2100   }
2101   else
2102   {
2103       ir->fepvals->n_lambda = 0;
2104   }
2105
2106   /* WALL PARAMETERS */
2107
2108   do_wall_params(ir,wall_atomtype,wall_density,opts);
2109
2110   /* ORIENTATION RESTRAINT PARAMETERS */
2111
2112   if (opts->bOrire && str_nelem(orirefitgrp,MAXPTR,NULL)!=1) {
2113       warning_error(wi,"ERROR: Need one orientation restraint fit group\n");
2114   }
2115
2116   /* DEFORMATION PARAMETERS */
2117
2118   clear_mat(ir->deform);
2119   for(i=0; i<6; i++)
2120   {
2121       dumdub[0][i] = 0;
2122   }
2123   m = sscanf(deform,"%lf %lf %lf %lf %lf %lf",
2124              &(dumdub[0][0]),&(dumdub[0][1]),&(dumdub[0][2]),
2125              &(dumdub[0][3]),&(dumdub[0][4]),&(dumdub[0][5]));
2126   for(i=0; i<3; i++)
2127   {
2128       ir->deform[i][i] = dumdub[0][i];
2129   }
2130   ir->deform[YY][XX] = dumdub[0][3];
2131   ir->deform[ZZ][XX] = dumdub[0][4];
2132   ir->deform[ZZ][YY] = dumdub[0][5];
2133   if (ir->epc != epcNO) {
2134     for(i=0; i<3; i++)
2135       for(j=0; j<=i; j++)
2136         if (ir->deform[i][j]!=0 && ir->compress[i][j]!=0) {
2137         warning_error(wi,"A box element has deform set and compressibility > 0");
2138         }
2139     for(i=0; i<3; i++)
2140       for(j=0; j<i; j++)
2141         if (ir->deform[i][j]!=0) {
2142           for(m=j; m<DIM; m++)
2143             if (ir->compress[m][j]!=0) {
2144               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.");
2145               warning(wi,warn_buf);
2146             }
2147         }
2148   }
2149
2150   sfree(dumstr[0]);
2151   sfree(dumstr[1]);
2152 }
2153
2154 static int search_QMstring(char *s,int ng,const char *gn[])
2155 {
2156   /* same as normal search_string, but this one searches QM strings */
2157   int i;
2158
2159   for(i=0; (i<ng); i++)
2160     if (gmx_strcasecmp(s,gn[i]) == 0)
2161       return i;
2162
2163   gmx_fatal(FARGS,"this QM method or basisset (%s) is not implemented\n!",s);
2164
2165   return -1;
2166
2167 } /* search_QMstring */
2168
2169
2170 int search_string(char *s,int ng,char *gn[])
2171 {
2172   int i;
2173
2174   for(i=0; (i<ng); i++)
2175   {
2176     if (gmx_strcasecmp(s,gn[i]) == 0)
2177     {
2178       return i;
2179     }
2180   }
2181
2182   gmx_fatal(FARGS,
2183             "Group %s referenced in the .mdp file was not found in the index file.\n"
2184             "Group names must match either [moleculetype] names or custom index group\n"
2185             "names, in which case you must supply an index file to the '-n' option\n"
2186             "of grompp.",
2187             s);
2188
2189   return -1;
2190 }
2191
2192 static gmx_bool do_numbering(int natoms,gmx_groups_t *groups,int ng,char *ptrs[],
2193                          t_blocka *block,char *gnames[],
2194                          int gtype,int restnm,
2195                          int grptp,gmx_bool bVerbose,
2196                          warninp_t wi)
2197 {
2198     unsigned short *cbuf;
2199     t_grps *grps=&(groups->grps[gtype]);
2200     int    i,j,gid,aj,ognr,ntot=0;
2201     const char *title;
2202     gmx_bool   bRest;
2203     char   warn_buf[STRLEN];
2204
2205     if (debug)
2206     {
2207         fprintf(debug,"Starting numbering %d groups of type %d\n",ng,gtype);
2208     }
2209
2210     title = gtypes[gtype];
2211
2212     snew(cbuf,natoms);
2213     /* Mark all id's as not set */
2214     for(i=0; (i<natoms); i++)
2215     {
2216         cbuf[i] = NOGID;
2217     }
2218
2219     snew(grps->nm_ind,ng+1); /* +1 for possible rest group */
2220     for(i=0; (i<ng); i++)
2221     {
2222         /* Lookup the group name in the block structure */
2223         gid = search_string(ptrs[i],block->nr,gnames);
2224         if ((grptp != egrptpONE) || (i == 0))
2225         {
2226             grps->nm_ind[grps->nr++]=gid;
2227         }
2228         if (debug)
2229         {
2230             fprintf(debug,"Found gid %d for group %s\n",gid,ptrs[i]);
2231         }
2232
2233         /* Now go over the atoms in the group */
2234         for(j=block->index[gid]; (j<block->index[gid+1]); j++)
2235         {
2236
2237             aj=block->a[j];
2238
2239             /* Range checking */
2240             if ((aj < 0) || (aj >= natoms))
2241             {
2242                 gmx_fatal(FARGS,"Invalid atom number %d in indexfile",aj);
2243             }
2244             /* Lookup up the old group number */
2245             ognr = cbuf[aj];
2246             if (ognr != NOGID)
2247             {
2248                 gmx_fatal(FARGS,"Atom %d in multiple %s groups (%d and %d)",
2249                           aj+1,title,ognr+1,i+1);
2250             }
2251             else
2252             {
2253                 /* Store the group number in buffer */
2254                 if (grptp == egrptpONE)
2255                 {
2256                     cbuf[aj] = 0;
2257                 }
2258                 else
2259                 {
2260                     cbuf[aj] = i;
2261                 }
2262                 ntot++;
2263             }
2264         }
2265     }
2266
2267     /* Now check whether we have done all atoms */
2268     bRest = FALSE;
2269     if (ntot != natoms)
2270     {
2271         if (grptp == egrptpALL)
2272         {
2273             gmx_fatal(FARGS,"%d atoms are not part of any of the %s groups",
2274                       natoms-ntot,title);
2275         }
2276         else if (grptp == egrptpPART)
2277         {
2278             sprintf(warn_buf,"%d atoms are not part of any of the %s groups",
2279                     natoms-ntot,title);
2280             warning_note(wi,warn_buf);
2281         }
2282         /* Assign all atoms currently unassigned to a rest group */
2283         for(j=0; (j<natoms); j++)
2284         {
2285             if (cbuf[j] == NOGID)
2286             {
2287                 cbuf[j] = grps->nr;
2288                 bRest = TRUE;
2289             }
2290         }
2291         if (grptp != egrptpPART)
2292         {
2293             if (bVerbose)
2294             {
2295                 fprintf(stderr,
2296                         "Making dummy/rest group for %s containing %d elements\n",
2297                         title,natoms-ntot);
2298             }
2299             /* Add group name "rest" */
2300             grps->nm_ind[grps->nr] = restnm;
2301
2302             /* Assign the rest name to all atoms not currently assigned to a group */
2303             for(j=0; (j<natoms); j++)
2304             {
2305                 if (cbuf[j] == NOGID)
2306                 {
2307                     cbuf[j] = grps->nr;
2308                 }
2309             }
2310             grps->nr++;
2311         }
2312     }
2313
2314     if (grps->nr == 1 && (ntot == 0 || ntot == natoms))
2315     {
2316         /* All atoms are part of one (or no) group, no index required */
2317         groups->ngrpnr[gtype] = 0;
2318         groups->grpnr[gtype]  = NULL;
2319     }
2320     else
2321     {
2322         groups->ngrpnr[gtype] = natoms;
2323         snew(groups->grpnr[gtype],natoms);
2324         for(j=0; (j<natoms); j++)
2325         {
2326             groups->grpnr[gtype][j] = cbuf[j];
2327         }
2328     }
2329
2330     sfree(cbuf);
2331
2332     return (bRest && grptp == egrptpPART);
2333 }
2334
2335 static void calc_nrdf(gmx_mtop_t *mtop,t_inputrec *ir,char **gnames)
2336 {
2337   t_grpopts *opts;
2338   gmx_groups_t *groups;
2339   t_pull  *pull;
2340   int     natoms,ai,aj,i,j,d,g,imin,jmin,nc;
2341   t_iatom *ia;
2342   int     *nrdf2,*na_vcm,na_tot;
2343   double  *nrdf_tc,*nrdf_vcm,nrdf_uc,n_sub=0;
2344   gmx_mtop_atomloop_all_t aloop;
2345   t_atom  *atom;
2346   int     mb,mol,ftype,as;
2347   gmx_molblock_t *molb;
2348   gmx_moltype_t *molt;
2349
2350   /* Calculate nrdf.
2351    * First calc 3xnr-atoms for each group
2352    * then subtract half a degree of freedom for each constraint
2353    *
2354    * Only atoms and nuclei contribute to the degrees of freedom...
2355    */
2356
2357   opts = &ir->opts;
2358
2359   groups = &mtop->groups;
2360   natoms = mtop->natoms;
2361
2362   /* Allocate one more for a possible rest group */
2363   /* We need to sum degrees of freedom into doubles,
2364    * since floats give too low nrdf's above 3 million atoms.
2365    */
2366   snew(nrdf_tc,groups->grps[egcTC].nr+1);
2367   snew(nrdf_vcm,groups->grps[egcVCM].nr+1);
2368   snew(na_vcm,groups->grps[egcVCM].nr+1);
2369
2370   for(i=0; i<groups->grps[egcTC].nr; i++)
2371     nrdf_tc[i] = 0;
2372   for(i=0; i<groups->grps[egcVCM].nr+1; i++)
2373     nrdf_vcm[i] = 0;
2374
2375   snew(nrdf2,natoms);
2376   aloop = gmx_mtop_atomloop_all_init(mtop);
2377   while (gmx_mtop_atomloop_all_next(aloop,&i,&atom)) {
2378     nrdf2[i] = 0;
2379     if (atom->ptype == eptAtom || atom->ptype == eptNucleus) {
2380       g = ggrpnr(groups,egcFREEZE,i);
2381       /* Double count nrdf for particle i */
2382       for(d=0; d<DIM; d++) {
2383         if (opts->nFreeze[g][d] == 0) {
2384           nrdf2[i] += 2;
2385         }
2386       }
2387       nrdf_tc [ggrpnr(groups,egcTC ,i)] += 0.5*nrdf2[i];
2388       nrdf_vcm[ggrpnr(groups,egcVCM,i)] += 0.5*nrdf2[i];
2389     }
2390   }
2391
2392   as = 0;
2393   for(mb=0; mb<mtop->nmolblock; mb++) {
2394     molb = &mtop->molblock[mb];
2395     molt = &mtop->moltype[molb->type];
2396     atom = molt->atoms.atom;
2397     for(mol=0; mol<molb->nmol; mol++) {
2398       for (ftype=F_CONSTR; ftype<=F_CONSTRNC; ftype++) {
2399         ia = molt->ilist[ftype].iatoms;
2400         for(i=0; i<molt->ilist[ftype].nr; ) {
2401           /* Subtract degrees of freedom for the constraints,
2402            * if the particles still have degrees of freedom left.
2403            * If one of the particles is a vsite or a shell, then all
2404            * constraint motion will go there, but since they do not
2405            * contribute to the constraints the degrees of freedom do not
2406            * change.
2407            */
2408           ai = as + ia[1];
2409           aj = as + ia[2];
2410           if (((atom[ia[1]].ptype == eptNucleus) ||
2411                (atom[ia[1]].ptype == eptAtom)) &&
2412               ((atom[ia[2]].ptype == eptNucleus) ||
2413                (atom[ia[2]].ptype == eptAtom))) {
2414             if (nrdf2[ai] > 0)
2415               jmin = 1;
2416             else
2417               jmin = 2;
2418             if (nrdf2[aj] > 0)
2419               imin = 1;
2420             else
2421               imin = 2;
2422             imin = min(imin,nrdf2[ai]);
2423             jmin = min(jmin,nrdf2[aj]);
2424             nrdf2[ai] -= imin;
2425             nrdf2[aj] -= jmin;
2426             nrdf_tc [ggrpnr(groups,egcTC ,ai)] -= 0.5*imin;
2427             nrdf_tc [ggrpnr(groups,egcTC ,aj)] -= 0.5*jmin;
2428             nrdf_vcm[ggrpnr(groups,egcVCM,ai)] -= 0.5*imin;
2429             nrdf_vcm[ggrpnr(groups,egcVCM,aj)] -= 0.5*jmin;
2430           }
2431           ia += interaction_function[ftype].nratoms+1;
2432           i  += interaction_function[ftype].nratoms+1;
2433         }
2434       }
2435       ia = molt->ilist[F_SETTLE].iatoms;
2436       for(i=0; i<molt->ilist[F_SETTLE].nr; ) {
2437         /* Subtract 1 dof from every atom in the SETTLE */
2438         for(j=0; j<3; j++) {
2439       ai = as + ia[1+j];
2440           imin = min(2,nrdf2[ai]);
2441           nrdf2[ai] -= imin;
2442           nrdf_tc [ggrpnr(groups,egcTC ,ai)] -= 0.5*imin;
2443           nrdf_vcm[ggrpnr(groups,egcVCM,ai)] -= 0.5*imin;
2444         }
2445         ia += 4;
2446         i  += 4;
2447       }
2448       as += molt->atoms.nr;
2449     }
2450   }
2451
2452   if (ir->ePull == epullCONSTRAINT) {
2453     /* Correct nrdf for the COM constraints.
2454      * We correct using the TC and VCM group of the first atom
2455      * in the reference and pull group. If atoms in one pull group
2456      * belong to different TC or VCM groups it is anyhow difficult
2457      * to determine the optimal nrdf assignment.
2458      */
2459     pull = ir->pull;
2460     if (pull->eGeom == epullgPOS) {
2461       nc = 0;
2462       for(i=0; i<DIM; i++) {
2463         if (pull->dim[i])
2464           nc++;
2465       }
2466     } else {
2467       nc = 1;
2468     }
2469     for(i=0; i<pull->ngrp; i++) {
2470       imin = 2*nc;
2471       if (pull->grp[0].nat > 0) {
2472         /* Subtract 1/2 dof from the reference group */
2473         ai = pull->grp[0].ind[0];
2474         if (nrdf_tc[ggrpnr(groups,egcTC,ai)] > 1) {
2475           nrdf_tc [ggrpnr(groups,egcTC ,ai)] -= 0.5;
2476           nrdf_vcm[ggrpnr(groups,egcVCM,ai)] -= 0.5;
2477           imin--;
2478         }
2479       }
2480       /* Subtract 1/2 dof from the pulled group */
2481       ai = pull->grp[1+i].ind[0];
2482       nrdf_tc [ggrpnr(groups,egcTC ,ai)] -= 0.5*imin;
2483       nrdf_vcm[ggrpnr(groups,egcVCM,ai)] -= 0.5*imin;
2484       if (nrdf_tc[ggrpnr(groups,egcTC,ai)] < 0)
2485         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)]]);
2486     }
2487   }
2488
2489   if (ir->nstcomm != 0) {
2490     /* Subtract 3 from the number of degrees of freedom in each vcm group
2491      * when com translation is removed and 6 when rotation is removed
2492      * as well.
2493      */
2494     switch (ir->comm_mode) {
2495     case ecmLINEAR:
2496       n_sub = ndof_com(ir);
2497       break;
2498     case ecmANGULAR:
2499       n_sub = 6;
2500       break;
2501     default:
2502       n_sub = 0;
2503       gmx_incons("Checking comm_mode");
2504     }
2505
2506     for(i=0; i<groups->grps[egcTC].nr; i++) {
2507       /* Count the number of atoms of TC group i for every VCM group */
2508       for(j=0; j<groups->grps[egcVCM].nr+1; j++)
2509         na_vcm[j] = 0;
2510       na_tot = 0;
2511       for(ai=0; ai<natoms; ai++)
2512         if (ggrpnr(groups,egcTC,ai) == i) {
2513           na_vcm[ggrpnr(groups,egcVCM,ai)]++;
2514           na_tot++;
2515         }
2516       /* Correct for VCM removal according to the fraction of each VCM
2517        * group present in this TC group.
2518        */
2519       nrdf_uc = nrdf_tc[i];
2520       if (debug) {
2521         fprintf(debug,"T-group[%d] nrdf_uc = %g, n_sub = %g\n",
2522                 i,nrdf_uc,n_sub);
2523       }
2524       nrdf_tc[i] = 0;
2525       for(j=0; j<groups->grps[egcVCM].nr+1; j++) {
2526         if (nrdf_vcm[j] > n_sub) {
2527           nrdf_tc[i] += nrdf_uc*((double)na_vcm[j]/(double)na_tot)*
2528             (nrdf_vcm[j] - n_sub)/nrdf_vcm[j];
2529         }
2530         if (debug) {
2531           fprintf(debug,"  nrdf_vcm[%d] = %g, nrdf = %g\n",
2532                   j,nrdf_vcm[j],nrdf_tc[i]);
2533         }
2534       }
2535     }
2536   }
2537   for(i=0; (i<groups->grps[egcTC].nr); i++) {
2538     opts->nrdf[i] = nrdf_tc[i];
2539     if (opts->nrdf[i] < 0)
2540       opts->nrdf[i] = 0;
2541     fprintf(stderr,
2542             "Number of degrees of freedom in T-Coupling group %s is %.2f\n",
2543             gnames[groups->grps[egcTC].nm_ind[i]],opts->nrdf[i]);
2544   }
2545
2546   sfree(nrdf2);
2547   sfree(nrdf_tc);
2548   sfree(nrdf_vcm);
2549   sfree(na_vcm);
2550 }
2551
2552 static void decode_cos(char *s,t_cosines *cosine,gmx_bool bTime)
2553 {
2554   char   *t;
2555   char   format[STRLEN],f1[STRLEN];
2556   double a,phi;
2557   int    i;
2558
2559   t=strdup(s);
2560   trim(t);
2561
2562   cosine->n=0;
2563   cosine->a=NULL;
2564   cosine->phi=NULL;
2565   if (strlen(t)) {
2566     sscanf(t,"%d",&(cosine->n));
2567     if (cosine->n <= 0) {
2568       cosine->n=0;
2569     } else {
2570       snew(cosine->a,cosine->n);
2571       snew(cosine->phi,cosine->n);
2572
2573       sprintf(format,"%%*d");
2574       for(i=0; (i<cosine->n); i++) {
2575         strcpy(f1,format);
2576         strcat(f1,"%lf%lf");
2577         if (sscanf(t,f1,&a,&phi) < 2)
2578           gmx_fatal(FARGS,"Invalid input for electric field shift: '%s'",t);
2579         cosine->a[i]=a;
2580         cosine->phi[i]=phi;
2581         strcat(format,"%*lf%*lf");
2582       }
2583     }
2584   }
2585   sfree(t);
2586 }
2587
2588 static gmx_bool do_egp_flag(t_inputrec *ir,gmx_groups_t *groups,
2589                         const char *option,const char *val,int flag)
2590 {
2591   /* The maximum number of energy group pairs would be MAXPTR*(MAXPTR+1)/2.
2592    * But since this is much larger than STRLEN, such a line can not be parsed.
2593    * The real maximum is the number of names that fit in a string: STRLEN/2.
2594    */
2595 #define EGP_MAX (STRLEN/2)
2596   int  nelem,i,j,k,nr;
2597   char *names[EGP_MAX];
2598   char ***gnames;
2599   gmx_bool bSet;
2600
2601   gnames = groups->grpname;
2602
2603   nelem = str_nelem(val,EGP_MAX,names);
2604   if (nelem % 2 != 0)
2605     gmx_fatal(FARGS,"The number of groups for %s is odd",option);
2606   nr = groups->grps[egcENER].nr;
2607   bSet = FALSE;
2608   for(i=0; i<nelem/2; i++) {
2609     j = 0;
2610     while ((j < nr) &&
2611            gmx_strcasecmp(names[2*i],*(gnames[groups->grps[egcENER].nm_ind[j]])))
2612       j++;
2613     if (j == nr)
2614       gmx_fatal(FARGS,"%s in %s is not an energy group\n",
2615                   names[2*i],option);
2616     k = 0;
2617     while ((k < nr) &&
2618            gmx_strcasecmp(names[2*i+1],*(gnames[groups->grps[egcENER].nm_ind[k]])))
2619       k++;
2620     if (k==nr)
2621       gmx_fatal(FARGS,"%s in %s is not an energy group\n",
2622               names[2*i+1],option);
2623     if ((j < nr) && (k < nr)) {
2624       ir->opts.egp_flags[nr*j+k] |= flag;
2625       ir->opts.egp_flags[nr*k+j] |= flag;
2626       bSet = TRUE;
2627     }
2628   }
2629
2630   return bSet;
2631 }
2632
2633 void do_index(const char* mdparin, const char *ndx,
2634               gmx_mtop_t *mtop,
2635               gmx_bool bVerbose,
2636               t_inputrec *ir,rvec *v,
2637               warninp_t wi)
2638 {
2639   t_blocka *grps;
2640   gmx_groups_t *groups;
2641   int     natoms;
2642   t_symtab *symtab;
2643   t_atoms atoms_all;
2644   char    warnbuf[STRLEN],**gnames;
2645   int     nr,ntcg,ntau_t,nref_t,nacc,nofg,nSA,nSA_points,nSA_time,nSA_temp;
2646   real    tau_min;
2647   int     nstcmin;
2648   int     nacg,nfreeze,nfrdim,nenergy,nvcm,nuser;
2649   char    *ptr1[MAXPTR],*ptr2[MAXPTR],*ptr3[MAXPTR];
2650   int     i,j,k,restnm;
2651   real    SAtime;
2652   gmx_bool    bExcl,bTable,bSetTCpar,bAnneal,bRest;
2653   int     nQMmethod,nQMbasis,nQMcharge,nQMmult,nbSH,nCASorb,nCASelec,
2654     nSAon,nSAoff,nSAsteps,nQMg,nbOPT,nbTS;
2655   char    warn_buf[STRLEN];
2656
2657   if (bVerbose)
2658     fprintf(stderr,"processing index file...\n");
2659   debug_gmx();
2660   if (ndx == NULL) {
2661     snew(grps,1);
2662     snew(grps->index,1);
2663     snew(gnames,1);
2664     atoms_all = gmx_mtop_global_atoms(mtop);
2665     analyse(&atoms_all,grps,&gnames,FALSE,TRUE);
2666     free_t_atoms(&atoms_all,FALSE);
2667   } else {
2668     grps = init_index(ndx,&gnames);
2669   }
2670
2671   groups = &mtop->groups;
2672   natoms = mtop->natoms;
2673   symtab = &mtop->symtab;
2674
2675   snew(groups->grpname,grps->nr+1);
2676
2677   for(i=0; (i<grps->nr); i++) {
2678     groups->grpname[i] = put_symtab(symtab,gnames[i]);
2679   }
2680   groups->grpname[i] = put_symtab(symtab,"rest");
2681   restnm=i;
2682   srenew(gnames,grps->nr+1);
2683   gnames[restnm] = *(groups->grpname[i]);
2684   groups->ngrpname = grps->nr+1;
2685
2686   set_warning_line(wi,mdparin,-1);
2687
2688   ntau_t = str_nelem(tau_t,MAXPTR,ptr1);
2689   nref_t = str_nelem(ref_t,MAXPTR,ptr2);
2690   ntcg   = str_nelem(tcgrps,MAXPTR,ptr3);
2691   if ((ntau_t != ntcg) || (nref_t != ntcg)) {
2692     gmx_fatal(FARGS,"Invalid T coupling input: %d groups, %d ref-t values and "
2693                 "%d tau-t values",ntcg,nref_t,ntau_t);
2694   }
2695
2696   bSetTCpar = (ir->etc || EI_SD(ir->eI) || ir->eI==eiBD || EI_TPI(ir->eI));
2697   do_numbering(natoms,groups,ntcg,ptr3,grps,gnames,egcTC,
2698                restnm,bSetTCpar ? egrptpALL : egrptpALL_GENREST,bVerbose,wi);
2699   nr = groups->grps[egcTC].nr;
2700   ir->opts.ngtc = nr;
2701   snew(ir->opts.nrdf,nr);
2702   snew(ir->opts.tau_t,nr);
2703   snew(ir->opts.ref_t,nr);
2704   if (ir->eI==eiBD && ir->bd_fric==0) {
2705     fprintf(stderr,"bd-fric=0, so tau-t will be used as the inverse friction constant(s)\n");
2706   }
2707
2708   if (bSetTCpar)
2709   {
2710       if (nr != nref_t)
2711       {
2712           gmx_fatal(FARGS,"Not enough ref-t and tau-t values!");
2713       }
2714
2715       tau_min = 1e20;
2716       for(i=0; (i<nr); i++)
2717       {
2718           ir->opts.tau_t[i] = strtod(ptr1[i],NULL);
2719           if ((ir->eI == eiBD || ir->eI == eiSD2) && ir->opts.tau_t[i] <= 0)
2720           {
2721               sprintf(warn_buf,"With integrator %s tau-t should be larger than 0",ei_names[ir->eI]);
2722               warning_error(wi,warn_buf);
2723           }
2724
2725           if (ir->etc != etcVRESCALE && ir->opts.tau_t[i] == 0)
2726           {
2727               warning_note(wi,"tau-t = -1 is the value to signal that a group should not have temperature coupling. Treating your use of tau-t = 0 as if you used -1.");
2728           }
2729
2730           if (ir->opts.tau_t[i] >= 0)
2731           {
2732               tau_min = min(tau_min,ir->opts.tau_t[i]);
2733           }
2734       }
2735       if (ir->etc != etcNO && ir->nsttcouple == -1)
2736       {
2737             ir->nsttcouple = ir_optimal_nsttcouple(ir);
2738       }
2739
2740       if (EI_VV(ir->eI))
2741       {
2742           if ((ir->etc==etcNOSEHOOVER) && (ir->epc==epcBERENDSEN)) {
2743               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");
2744           }
2745           if ((ir->epc==epcMTTK) && (ir->etc>etcNO))
2746           {
2747               if (ir->nstpcouple != ir->nsttcouple)
2748               {
2749                   int mincouple = min(ir->nstpcouple,ir->nsttcouple);
2750                   ir->nstpcouple = ir->nsttcouple = mincouple;
2751                   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);
2752                   warning_note(wi,warn_buf);
2753               }
2754           }
2755       }
2756       /* velocity verlet with averaged kinetic energy KE = 0.5*(v(t+1/2) - v(t-1/2)) is implemented
2757          primarily for testing purposes, and does not work with temperature coupling other than 1 */
2758
2759       if (ETC_ANDERSEN(ir->etc)) {
2760           if (ir->nsttcouple != 1) {
2761               ir->nsttcouple = 1;
2762               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");
2763               warning_note(wi,warn_buf);
2764           }
2765       }
2766       nstcmin = tcouple_min_integration_steps(ir->etc);
2767       if (nstcmin > 1)
2768       {
2769           if (tau_min/(ir->delta_t*ir->nsttcouple) < nstcmin)
2770           {
2771               sprintf(warn_buf,"For proper integration of the %s thermostat, tau-t (%g) should be at least %d times larger than nsttcouple*dt (%g)",
2772                       ETCOUPLTYPE(ir->etc),
2773                       tau_min,nstcmin,
2774                       ir->nsttcouple*ir->delta_t);
2775               warning(wi,warn_buf);
2776           }
2777       }
2778       for(i=0; (i<nr); i++)
2779       {
2780           ir->opts.ref_t[i] = strtod(ptr2[i],NULL);
2781           if (ir->opts.ref_t[i] < 0)
2782           {
2783               gmx_fatal(FARGS,"ref-t for group %d negative",i);
2784           }
2785       }
2786       /* set the lambda mc temperature to the md integrator temperature (which should be defined
2787          if we are in this conditional) if mc_temp is negative */
2788       if (ir->expandedvals->mc_temp < 0)
2789       {
2790           ir->expandedvals->mc_temp = ir->opts.ref_t[0];  /*for now, set to the first reft */
2791       }
2792   }
2793
2794   /* Simulated annealing for each group. There are nr groups */
2795   nSA = str_nelem(anneal,MAXPTR,ptr1);
2796   if (nSA == 1 && (ptr1[0][0]=='n' || ptr1[0][0]=='N'))
2797      nSA = 0;
2798   if(nSA>0 && nSA != nr)
2799     gmx_fatal(FARGS,"Not enough annealing values: %d (for %d groups)\n",nSA,nr);
2800   else {
2801     snew(ir->opts.annealing,nr);
2802     snew(ir->opts.anneal_npoints,nr);
2803     snew(ir->opts.anneal_time,nr);
2804     snew(ir->opts.anneal_temp,nr);
2805     for(i=0;i<nr;i++) {
2806       ir->opts.annealing[i]=eannNO;
2807       ir->opts.anneal_npoints[i]=0;
2808       ir->opts.anneal_time[i]=NULL;
2809       ir->opts.anneal_temp[i]=NULL;
2810     }
2811     if (nSA > 0) {
2812       bAnneal=FALSE;
2813       for(i=0;i<nr;i++) {
2814         if(ptr1[i][0]=='n' || ptr1[i][0]=='N') {
2815           ir->opts.annealing[i]=eannNO;
2816         } else if(ptr1[i][0]=='s'|| ptr1[i][0]=='S') {
2817           ir->opts.annealing[i]=eannSINGLE;
2818           bAnneal=TRUE;
2819         } else if(ptr1[i][0]=='p'|| ptr1[i][0]=='P') {
2820           ir->opts.annealing[i]=eannPERIODIC;
2821           bAnneal=TRUE;
2822         }
2823       }
2824       if(bAnneal) {
2825         /* Read the other fields too */
2826         nSA_points = str_nelem(anneal_npoints,MAXPTR,ptr1);
2827         if(nSA_points!=nSA)
2828           gmx_fatal(FARGS,"Found %d annealing-npoints values for %d groups\n",nSA_points,nSA);
2829         for(k=0,i=0;i<nr;i++) {
2830           ir->opts.anneal_npoints[i]=strtol(ptr1[i],NULL,10);
2831           if(ir->opts.anneal_npoints[i]==1)
2832             gmx_fatal(FARGS,"Please specify at least a start and an end point for annealing\n");
2833           snew(ir->opts.anneal_time[i],ir->opts.anneal_npoints[i]);
2834           snew(ir->opts.anneal_temp[i],ir->opts.anneal_npoints[i]);
2835           k += ir->opts.anneal_npoints[i];
2836         }
2837
2838         nSA_time = str_nelem(anneal_time,MAXPTR,ptr1);
2839         if(nSA_time!=k)
2840           gmx_fatal(FARGS,"Found %d annealing-time values, wanter %d\n",nSA_time,k);
2841         nSA_temp = str_nelem(anneal_temp,MAXPTR,ptr2);
2842         if(nSA_temp!=k)
2843           gmx_fatal(FARGS,"Found %d annealing-temp values, wanted %d\n",nSA_temp,k);
2844
2845         for(i=0,k=0;i<nr;i++) {
2846
2847           for(j=0;j<ir->opts.anneal_npoints[i];j++) {
2848             ir->opts.anneal_time[i][j]=strtod(ptr1[k],NULL);
2849             ir->opts.anneal_temp[i][j]=strtod(ptr2[k],NULL);
2850             if(j==0) {
2851               if(ir->opts.anneal_time[i][0] > (ir->init_t+GMX_REAL_EPS))
2852                 gmx_fatal(FARGS,"First time point for annealing > init_t.\n");
2853             } else {
2854               /* j>0 */
2855               if(ir->opts.anneal_time[i][j]<ir->opts.anneal_time[i][j-1])
2856                 gmx_fatal(FARGS,"Annealing timepoints out of order: t=%f comes after t=%f\n",
2857                             ir->opts.anneal_time[i][j],ir->opts.anneal_time[i][j-1]);
2858             }
2859             if(ir->opts.anneal_temp[i][j]<0)
2860               gmx_fatal(FARGS,"Found negative temperature in annealing: %f\n",ir->opts.anneal_temp[i][j]);
2861             k++;
2862           }
2863         }
2864         /* Print out some summary information, to make sure we got it right */
2865         for(i=0,k=0;i<nr;i++) {
2866           if(ir->opts.annealing[i]!=eannNO) {
2867             j = groups->grps[egcTC].nm_ind[i];
2868             fprintf(stderr,"Simulated annealing for group %s: %s, %d timepoints\n",
2869                     *(groups->grpname[j]),eann_names[ir->opts.annealing[i]],
2870                     ir->opts.anneal_npoints[i]);
2871             fprintf(stderr,"Time (ps)   Temperature (K)\n");
2872             /* All terms except the last one */
2873             for(j=0;j<(ir->opts.anneal_npoints[i]-1);j++)
2874                 fprintf(stderr,"%9.1f      %5.1f\n",ir->opts.anneal_time[i][j],ir->opts.anneal_temp[i][j]);
2875
2876             /* Finally the last one */
2877             j = ir->opts.anneal_npoints[i]-1;
2878             if(ir->opts.annealing[i]==eannSINGLE)
2879               fprintf(stderr,"%9.1f-     %5.1f\n",ir->opts.anneal_time[i][j],ir->opts.anneal_temp[i][j]);
2880             else {
2881               fprintf(stderr,"%9.1f      %5.1f\n",ir->opts.anneal_time[i][j],ir->opts.anneal_temp[i][j]);
2882               if(fabs(ir->opts.anneal_temp[i][j]-ir->opts.anneal_temp[i][0])>GMX_REAL_EPS)
2883                 warning_note(wi,"There is a temperature jump when your annealing loops back.\n");
2884             }
2885           }
2886         }
2887       }
2888     }
2889   }
2890
2891   if (ir->ePull != epullNO) {
2892     make_pull_groups(ir->pull,pull_grp,grps,gnames);
2893   }
2894
2895   if (ir->bRot) {
2896     make_rotation_groups(ir->rot,rot_grp,grps,gnames);
2897   }
2898
2899   nacc = str_nelem(acc,MAXPTR,ptr1);
2900   nacg = str_nelem(accgrps,MAXPTR,ptr2);
2901   if (nacg*DIM != nacc)
2902     gmx_fatal(FARGS,"Invalid Acceleration input: %d groups and %d acc. values",
2903                 nacg,nacc);
2904   do_numbering(natoms,groups,nacg,ptr2,grps,gnames,egcACC,
2905                restnm,egrptpALL_GENREST,bVerbose,wi);
2906   nr = groups->grps[egcACC].nr;
2907   snew(ir->opts.acc,nr);
2908   ir->opts.ngacc=nr;
2909
2910   for(i=k=0; (i<nacg); i++)
2911     for(j=0; (j<DIM); j++,k++)
2912       ir->opts.acc[i][j]=strtod(ptr1[k],NULL);
2913   for( ;(i<nr); i++)
2914     for(j=0; (j<DIM); j++)
2915       ir->opts.acc[i][j]=0;
2916
2917   nfrdim  = str_nelem(frdim,MAXPTR,ptr1);
2918   nfreeze = str_nelem(freeze,MAXPTR,ptr2);
2919   if (nfrdim != DIM*nfreeze)
2920     gmx_fatal(FARGS,"Invalid Freezing input: %d groups and %d freeze values",
2921                 nfreeze,nfrdim);
2922   do_numbering(natoms,groups,nfreeze,ptr2,grps,gnames,egcFREEZE,
2923                restnm,egrptpALL_GENREST,bVerbose,wi);
2924   nr = groups->grps[egcFREEZE].nr;
2925   ir->opts.ngfrz=nr;
2926   snew(ir->opts.nFreeze,nr);
2927   for(i=k=0; (i<nfreeze); i++)
2928     for(j=0; (j<DIM); j++,k++) {
2929       ir->opts.nFreeze[i][j]=(gmx_strncasecmp(ptr1[k],"Y",1)==0);
2930       if (!ir->opts.nFreeze[i][j]) {
2931         if (gmx_strncasecmp(ptr1[k],"N",1) != 0) {
2932           sprintf(warnbuf,"Please use Y(ES) or N(O) for freezedim only "
2933                   "(not %s)", ptr1[k]);
2934           warning(wi,warn_buf);
2935         }
2936       }
2937     }
2938   for( ; (i<nr); i++)
2939     for(j=0; (j<DIM); j++)
2940       ir->opts.nFreeze[i][j]=0;
2941
2942   nenergy=str_nelem(energy,MAXPTR,ptr1);
2943   do_numbering(natoms,groups,nenergy,ptr1,grps,gnames,egcENER,
2944                restnm,egrptpALL_GENREST,bVerbose,wi);
2945   add_wall_energrps(groups,ir->nwall,symtab);
2946   ir->opts.ngener = groups->grps[egcENER].nr;
2947   nvcm=str_nelem(vcm,MAXPTR,ptr1);
2948   bRest =
2949     do_numbering(natoms,groups,nvcm,ptr1,grps,gnames,egcVCM,
2950                  restnm,nvcm==0 ? egrptpALL_GENREST : egrptpPART,bVerbose,wi);
2951   if (bRest) {
2952     warning(wi,"Some atoms are not part of any center of mass motion removal group.\n"
2953             "This may lead to artifacts.\n"
2954             "In most cases one should use one group for the whole system.");
2955   }
2956
2957   /* Now we have filled the freeze struct, so we can calculate NRDF */
2958   calc_nrdf(mtop,ir,gnames);
2959
2960   if (v && NULL) {
2961     real fac,ntot=0;
2962
2963     /* Must check per group! */
2964     for(i=0; (i<ir->opts.ngtc); i++)
2965       ntot += ir->opts.nrdf[i];
2966     if (ntot != (DIM*natoms)) {
2967       fac = sqrt(ntot/(DIM*natoms));
2968       if (bVerbose)
2969         fprintf(stderr,"Scaling velocities by a factor of %.3f to account for constraints\n"
2970                 "and removal of center of mass motion\n",fac);
2971       for(i=0; (i<natoms); i++)
2972         svmul(fac,v[i],v[i]);
2973     }
2974   }
2975
2976   nuser=str_nelem(user1,MAXPTR,ptr1);
2977   do_numbering(natoms,groups,nuser,ptr1,grps,gnames,egcUser1,
2978                restnm,egrptpALL_GENREST,bVerbose,wi);
2979   nuser=str_nelem(user2,MAXPTR,ptr1);
2980   do_numbering(natoms,groups,nuser,ptr1,grps,gnames,egcUser2,
2981                restnm,egrptpALL_GENREST,bVerbose,wi);
2982   nuser=str_nelem(xtc_grps,MAXPTR,ptr1);
2983   do_numbering(natoms,groups,nuser,ptr1,grps,gnames,egcXTC,
2984                restnm,egrptpONE,bVerbose,wi);
2985   nofg = str_nelem(orirefitgrp,MAXPTR,ptr1);
2986   do_numbering(natoms,groups,nofg,ptr1,grps,gnames,egcORFIT,
2987                restnm,egrptpALL_GENREST,bVerbose,wi);
2988
2989   /* QMMM input processing */
2990   nQMg          = str_nelem(QMMM,MAXPTR,ptr1);
2991   nQMmethod     = str_nelem(QMmethod,MAXPTR,ptr2);
2992   nQMbasis      = str_nelem(QMbasis,MAXPTR,ptr3);
2993   if((nQMmethod != nQMg)||(nQMbasis != nQMg)){
2994     gmx_fatal(FARGS,"Invalid QMMM input: %d groups %d basissets"
2995               " and %d methods\n",nQMg,nQMbasis,nQMmethod);
2996   }
2997   /* group rest, if any, is always MM! */
2998   do_numbering(natoms,groups,nQMg,ptr1,grps,gnames,egcQMMM,
2999                restnm,egrptpALL_GENREST,bVerbose,wi);
3000   nr = nQMg; /*atoms->grps[egcQMMM].nr;*/
3001   ir->opts.ngQM = nQMg;
3002   snew(ir->opts.QMmethod,nr);
3003   snew(ir->opts.QMbasis,nr);
3004   for(i=0;i<nr;i++){
3005     /* input consists of strings: RHF CASSCF PM3 .. These need to be
3006      * converted to the corresponding enum in names.c
3007      */
3008     ir->opts.QMmethod[i] = search_QMstring(ptr2[i],eQMmethodNR,
3009                                            eQMmethod_names);
3010     ir->opts.QMbasis[i]  = search_QMstring(ptr3[i],eQMbasisNR,
3011                                            eQMbasis_names);
3012
3013   }
3014   nQMmult   = str_nelem(QMmult,MAXPTR,ptr1);
3015   nQMcharge = str_nelem(QMcharge,MAXPTR,ptr2);
3016   nbSH      = str_nelem(bSH,MAXPTR,ptr3);
3017   snew(ir->opts.QMmult,nr);
3018   snew(ir->opts.QMcharge,nr);
3019   snew(ir->opts.bSH,nr);
3020
3021   for(i=0;i<nr;i++){
3022     ir->opts.QMmult[i]   = strtol(ptr1[i],NULL,10);
3023     ir->opts.QMcharge[i] = strtol(ptr2[i],NULL,10);
3024     ir->opts.bSH[i]      = (gmx_strncasecmp(ptr3[i],"Y",1)==0);
3025   }
3026
3027   nCASelec  = str_nelem(CASelectrons,MAXPTR,ptr1);
3028   nCASorb   = str_nelem(CASorbitals,MAXPTR,ptr2);
3029   snew(ir->opts.CASelectrons,nr);
3030   snew(ir->opts.CASorbitals,nr);
3031   for(i=0;i<nr;i++){
3032     ir->opts.CASelectrons[i]= strtol(ptr1[i],NULL,10);
3033     ir->opts.CASorbitals[i] = strtol(ptr2[i],NULL,10);
3034   }
3035   /* special optimization options */
3036
3037   nbOPT = str_nelem(bOPT,MAXPTR,ptr1);
3038   nbTS = str_nelem(bTS,MAXPTR,ptr2);
3039   snew(ir->opts.bOPT,nr);
3040   snew(ir->opts.bTS,nr);
3041   for(i=0;i<nr;i++){
3042     ir->opts.bOPT[i] = (gmx_strncasecmp(ptr1[i],"Y",1)==0);
3043     ir->opts.bTS[i]  = (gmx_strncasecmp(ptr2[i],"Y",1)==0);
3044   }
3045   nSAon     = str_nelem(SAon,MAXPTR,ptr1);
3046   nSAoff    = str_nelem(SAoff,MAXPTR,ptr2);
3047   nSAsteps  = str_nelem(SAsteps,MAXPTR,ptr3);
3048   snew(ir->opts.SAon,nr);
3049   snew(ir->opts.SAoff,nr);
3050   snew(ir->opts.SAsteps,nr);
3051
3052   for(i=0;i<nr;i++){
3053     ir->opts.SAon[i]    = strtod(ptr1[i],NULL);
3054     ir->opts.SAoff[i]   = strtod(ptr2[i],NULL);
3055     ir->opts.SAsteps[i] = strtol(ptr3[i],NULL,10);
3056   }
3057   /* end of QMMM input */
3058
3059   if (bVerbose)
3060     for(i=0; (i<egcNR); i++) {
3061       fprintf(stderr,"%-16s has %d element(s):",gtypes[i],groups->grps[i].nr);
3062       for(j=0; (j<groups->grps[i].nr); j++)
3063         fprintf(stderr," %s",*(groups->grpname[groups->grps[i].nm_ind[j]]));
3064       fprintf(stderr,"\n");
3065     }
3066
3067   nr = groups->grps[egcENER].nr;
3068   snew(ir->opts.egp_flags,nr*nr);
3069
3070   bExcl = do_egp_flag(ir,groups,"energygrp-excl",egpexcl,EGP_EXCL);
3071     if (bExcl && ir->cutoff_scheme == ecutsVERLET)
3072     {
3073         warning_error(wi,"Energy group exclusions are not (yet) implemented for the Verlet scheme");
3074     }
3075   if (bExcl && EEL_FULL(ir->coulombtype))
3076     warning(wi,"Can not exclude the lattice Coulomb energy between energy groups");
3077
3078   bTable = do_egp_flag(ir,groups,"energygrp-table",egptable,EGP_TABLE);
3079   if (bTable && !(ir->vdwtype == evdwUSER) &&
3080       !(ir->coulombtype == eelUSER) && !(ir->coulombtype == eelPMEUSER) &&
3081       !(ir->coulombtype == eelPMEUSERSWITCH))
3082     gmx_fatal(FARGS,"Can only have energy group pair tables in combination with user tables for VdW and/or Coulomb");
3083
3084   decode_cos(efield_x,&(ir->ex[XX]),FALSE);
3085   decode_cos(efield_xt,&(ir->et[XX]),TRUE);
3086   decode_cos(efield_y,&(ir->ex[YY]),FALSE);
3087   decode_cos(efield_yt,&(ir->et[YY]),TRUE);
3088   decode_cos(efield_z,&(ir->ex[ZZ]),FALSE);
3089   decode_cos(efield_zt,&(ir->et[ZZ]),TRUE);
3090
3091   if (ir->bAdress)
3092     do_adress_index(ir->adress,groups,gnames,&(ir->opts),wi);
3093
3094   for(i=0; (i<grps->nr); i++)
3095     sfree(gnames[i]);
3096   sfree(gnames);
3097   done_blocka(grps);
3098   sfree(grps);
3099
3100 }
3101
3102
3103
3104 static void check_disre(gmx_mtop_t *mtop)
3105 {
3106   gmx_ffparams_t *ffparams;
3107   t_functype *functype;
3108   t_iparams  *ip;
3109   int i,ndouble,ftype;
3110   int label,old_label;
3111
3112   if (gmx_mtop_ftype_count(mtop,F_DISRES) > 0) {
3113     ffparams  = &mtop->ffparams;
3114     functype  = ffparams->functype;
3115     ip        = ffparams->iparams;
3116     ndouble   = 0;
3117     old_label = -1;
3118     for(i=0; i<ffparams->ntypes; i++) {
3119       ftype = functype[i];
3120       if (ftype == F_DISRES) {
3121         label = ip[i].disres.label;
3122         if (label == old_label) {
3123           fprintf(stderr,"Distance restraint index %d occurs twice\n",label);
3124           ndouble++;
3125         }
3126         old_label = label;
3127       }
3128     }
3129     if (ndouble>0)
3130       gmx_fatal(FARGS,"Found %d double distance restraint indices,\n"
3131                 "probably the parameters for multiple pairs in one restraint "
3132                 "are not identical\n",ndouble);
3133   }
3134 }
3135
3136 static gmx_bool absolute_reference(t_inputrec *ir,gmx_mtop_t *sys,
3137                                    gmx_bool posres_only,
3138                                    ivec AbsRef)
3139 {
3140     int d,g,i;
3141     gmx_mtop_ilistloop_t iloop;
3142     t_ilist *ilist;
3143     int nmol;
3144     t_iparams *pr;
3145
3146     clear_ivec(AbsRef);
3147
3148     if (!posres_only)
3149     {
3150         /* Check the COM */
3151         for(d=0; d<DIM; d++)
3152         {
3153             AbsRef[d] = (d < ndof_com(ir) ? 0 : 1);
3154         }
3155         /* Check for freeze groups */
3156         for(g=0; g<ir->opts.ngfrz; g++)
3157         {
3158             for(d=0; d<DIM; d++)
3159             {
3160                 if (ir->opts.nFreeze[g][d] != 0)
3161                 {
3162                     AbsRef[d] = 1;
3163                 }
3164             }
3165         }
3166     }
3167
3168     /* Check for position restraints */
3169     iloop = gmx_mtop_ilistloop_init(sys);
3170     while (gmx_mtop_ilistloop_next(iloop,&ilist,&nmol))
3171     {
3172         if (nmol > 0 &&
3173             (AbsRef[XX] == 0 || AbsRef[YY] == 0 || AbsRef[ZZ] == 0))
3174         {
3175             for(i=0; i<ilist[F_POSRES].nr; i+=2)
3176             {
3177                 pr = &sys->ffparams.iparams[ilist[F_POSRES].iatoms[i]];
3178                 for(d=0; d<DIM; d++)
3179                 {
3180                     if (pr->posres.fcA[d] != 0)
3181                     {
3182                         AbsRef[d] = 1;
3183                     }
3184                 }
3185             }
3186         }
3187     }
3188
3189     return (AbsRef[XX] != 0 && AbsRef[YY] != 0 && AbsRef[ZZ] != 0);
3190 }
3191
3192 void triple_check(const char *mdparin,t_inputrec *ir,gmx_mtop_t *sys,
3193                   warninp_t wi)
3194 {
3195   char err_buf[256];
3196   int  i,m,g,nmol,npct;
3197   gmx_bool bCharge,bAcc;
3198   real gdt_max,*mgrp,mt;
3199   rvec acc;
3200   gmx_mtop_atomloop_block_t aloopb;
3201   gmx_mtop_atomloop_all_t aloop;
3202   t_atom *atom;
3203   ivec AbsRef;
3204   char warn_buf[STRLEN];
3205
3206   set_warning_line(wi,mdparin,-1);
3207
3208   if (EI_DYNAMICS(ir->eI) && !EI_SD(ir->eI) && ir->eI != eiBD &&
3209       ir->comm_mode == ecmNO &&
3210       !(absolute_reference(ir,sys,FALSE,AbsRef) || ir->nsteps <= 10)) {
3211     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");
3212   }
3213
3214     /* Check for pressure coupling with absolute position restraints */
3215     if (ir->epc != epcNO && ir->refcoord_scaling == erscNO)
3216     {
3217         absolute_reference(ir,sys,TRUE,AbsRef);
3218         {
3219             for(m=0; m<DIM; m++)
3220             {
3221                 if (AbsRef[m] && norm2(ir->compress[m]) > 0)
3222                 {
3223                     warning(wi,"You are using pressure coupling with absolute position restraints, this will give artifacts. Use the refcoord_scaling option.");
3224                     break;
3225                 }
3226             }
3227         }
3228     }
3229
3230   bCharge = FALSE;
3231   aloopb = gmx_mtop_atomloop_block_init(sys);
3232   while (gmx_mtop_atomloop_block_next(aloopb,&atom,&nmol)) {
3233     if (atom->q != 0 || atom->qB != 0) {
3234       bCharge = TRUE;
3235     }
3236   }
3237
3238   if (!bCharge) {
3239     if (EEL_FULL(ir->coulombtype)) {
3240       sprintf(err_buf,
3241               "You are using full electrostatics treatment %s for a system without charges.\n"
3242               "This costs a lot of performance for just processing zeros, consider using %s instead.\n",
3243               EELTYPE(ir->coulombtype),EELTYPE(eelCUT));
3244       warning(wi,err_buf);
3245     }
3246   } else {
3247     if (ir->coulombtype == eelCUT && ir->rcoulomb > 0 && !ir->implicit_solvent) {
3248       sprintf(err_buf,
3249               "You are using a plain Coulomb cut-off, which might produce artifacts.\n"
3250               "You might want to consider using %s electrostatics.\n",
3251               EELTYPE(eelPME));
3252       warning_note(wi,err_buf);
3253     }
3254   }
3255
3256   /* Generalized reaction field */
3257   if (ir->opts.ngtc == 0) {
3258     sprintf(err_buf,"No temperature coupling while using coulombtype %s",
3259             eel_names[eelGRF]);
3260     CHECK(ir->coulombtype == eelGRF);
3261   }
3262   else {
3263     sprintf(err_buf,"When using coulombtype = %s"
3264             " ref-t for temperature coupling should be > 0",
3265             eel_names[eelGRF]);
3266     CHECK((ir->coulombtype == eelGRF) && (ir->opts.ref_t[0] <= 0));
3267   }
3268
3269     if (ir->eI == eiSD1 &&
3270         (gmx_mtop_ftype_count(sys,F_CONSTR) > 0 ||
3271          gmx_mtop_ftype_count(sys,F_SETTLE) > 0))
3272     {
3273         sprintf(warn_buf,"With constraints integrator %s is less accurate, consider using %s instead",ei_names[ir->eI],ei_names[eiSD2]);
3274         warning_note(wi,warn_buf);
3275     }
3276
3277   bAcc = FALSE;
3278   for(i=0; (i<sys->groups.grps[egcACC].nr); i++) {
3279     for(m=0; (m<DIM); m++) {
3280       if (fabs(ir->opts.acc[i][m]) > 1e-6) {
3281         bAcc = TRUE;
3282       }
3283     }
3284   }
3285   if (bAcc) {
3286     clear_rvec(acc);
3287     snew(mgrp,sys->groups.grps[egcACC].nr);
3288     aloop = gmx_mtop_atomloop_all_init(sys);
3289     while (gmx_mtop_atomloop_all_next(aloop,&i,&atom)) {
3290       mgrp[ggrpnr(&sys->groups,egcACC,i)] += atom->m;
3291     }
3292     mt = 0.0;
3293     for(i=0; (i<sys->groups.grps[egcACC].nr); i++) {
3294       for(m=0; (m<DIM); m++)
3295         acc[m] += ir->opts.acc[i][m]*mgrp[i];
3296       mt += mgrp[i];
3297     }
3298     for(m=0; (m<DIM); m++) {
3299       if (fabs(acc[m]) > 1e-6) {
3300         const char *dim[DIM] = { "X", "Y", "Z" };
3301         fprintf(stderr,
3302                 "Net Acceleration in %s direction, will %s be corrected\n",
3303                 dim[m],ir->nstcomm != 0 ? "" : "not");
3304         if (ir->nstcomm != 0 && m < ndof_com(ir)) {
3305           acc[m] /= mt;
3306           for (i=0; (i<sys->groups.grps[egcACC].nr); i++)
3307             ir->opts.acc[i][m] -= acc[m];
3308         }
3309       }
3310     }
3311     sfree(mgrp);
3312   }
3313
3314   if (ir->efep != efepNO && ir->fepvals->sc_alpha != 0 &&
3315       !gmx_within_tol(sys->ffparams.reppow,12.0,10*GMX_DOUBLE_EPS)) {
3316     gmx_fatal(FARGS,"Soft-core interactions are only supported with VdW repulsion power 12");
3317   }
3318
3319   if (ir->ePull != epullNO) {
3320     if (ir->pull->grp[0].nat == 0) {
3321         absolute_reference(ir,sys,FALSE,AbsRef);
3322       for(m=0; m<DIM; m++) {
3323         if (ir->pull->dim[m] && !AbsRef[m]) {
3324           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.");
3325           break;
3326         }
3327       }
3328     }
3329
3330     if (ir->pull->eGeom == epullgDIRPBC) {
3331       for(i=0; i<3; i++) {
3332         for(m=0; m<=i; m++) {
3333           if ((ir->epc != epcNO && ir->compress[i][m] != 0) ||
3334               ir->deform[i][m] != 0) {
3335             for(g=1; g<ir->pull->ngrp; g++) {
3336               if (ir->pull->grp[g].vec[m] != 0) {
3337                 gmx_fatal(FARGS,"Can not have dynamic box while using pull geometry '%s' (dim %c)",EPULLGEOM(ir->pull->eGeom),'x'+m);
3338               }
3339             }
3340           }
3341         }
3342       }
3343     }
3344   }
3345
3346   check_disre(sys);
3347 }
3348
3349 void double_check(t_inputrec *ir,matrix box,gmx_bool bConstr,warninp_t wi)
3350 {
3351   real min_size;
3352   gmx_bool bTWIN;
3353   char warn_buf[STRLEN];
3354   const char *ptr;
3355
3356   ptr = check_box(ir->ePBC,box);
3357   if (ptr) {
3358       warning_error(wi,ptr);
3359   }
3360
3361   if (bConstr && ir->eConstrAlg == econtSHAKE) {
3362     if (ir->shake_tol <= 0.0) {
3363       sprintf(warn_buf,"ERROR: shake-tol must be > 0 instead of %g\n",
3364               ir->shake_tol);
3365       warning_error(wi,warn_buf);
3366     }
3367
3368     if (IR_TWINRANGE(*ir) && ir->nstlist > 1) {
3369       sprintf(warn_buf,"With twin-range cut-off's and SHAKE the virial and the pressure are incorrect.");
3370       if (ir->epc == epcNO) {
3371         warning(wi,warn_buf);
3372       } else {
3373           warning_error(wi,warn_buf);
3374       }
3375     }
3376   }
3377
3378   if( (ir->eConstrAlg == econtLINCS) && bConstr) {
3379     /* If we have Lincs constraints: */
3380     if(ir->eI==eiMD && ir->etc==etcNO &&
3381        ir->eConstrAlg==econtLINCS && ir->nLincsIter==1) {
3382       sprintf(warn_buf,"For energy conservation with LINCS, lincs_iter should be 2 or larger.\n");
3383       warning_note(wi,warn_buf);
3384     }
3385
3386     if ((ir->eI == eiCG || ir->eI == eiLBFGS) && (ir->nProjOrder<8)) {
3387       sprintf(warn_buf,"For accurate %s with LINCS constraints, lincs-order should be 8 or more.",ei_names[ir->eI]);
3388       warning_note(wi,warn_buf);
3389     }
3390     if (ir->epc==epcMTTK) {
3391         warning_error(wi,"MTTK not compatible with lincs -- use shake instead.");
3392     }
3393   }
3394
3395   if (ir->LincsWarnAngle > 90.0) {
3396     sprintf(warn_buf,"lincs-warnangle can not be larger than 90 degrees, setting it to 90.\n");
3397     warning(wi,warn_buf);
3398     ir->LincsWarnAngle = 90.0;
3399   }
3400
3401   if (ir->ePBC != epbcNONE) {
3402     if (ir->nstlist == 0) {
3403       warning(wi,"With nstlist=0 atoms are only put into the box at step 0, therefore drifting atoms might cause the simulation to crash.");
3404     }
3405     bTWIN = (ir->rlistlong > ir->rlist);
3406     if (ir->ns_type == ensGRID) {
3407       if (sqr(ir->rlistlong) >= max_cutoff2(ir->ePBC,box)) {
3408           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",
3409                 bTWIN ? (ir->rcoulomb==ir->rlistlong ? "rcoulomb" : "rvdw"):"rlist");
3410           warning_error(wi,warn_buf);
3411       }
3412     } else {
3413       min_size = min(box[XX][XX],min(box[YY][YY],box[ZZ][ZZ]));
3414       if (2*ir->rlistlong >= min_size) {
3415           sprintf(warn_buf,"ERROR: One of the box lengths is smaller than twice the cut-off length. Increase the box size or decrease rlist.");
3416           warning_error(wi,warn_buf);
3417         if (TRICLINIC(box))
3418           fprintf(stderr,"Grid search might allow larger cut-off's than simple search with triclinic boxes.");
3419       }
3420     }
3421   }
3422 }
3423
3424 void check_chargegroup_radii(const gmx_mtop_t *mtop,const t_inputrec *ir,
3425                              rvec *x,
3426                              warninp_t wi)
3427 {
3428     real rvdw1,rvdw2,rcoul1,rcoul2;
3429     char warn_buf[STRLEN];
3430
3431     calc_chargegroup_radii(mtop,x,&rvdw1,&rvdw2,&rcoul1,&rcoul2);
3432
3433     if (rvdw1 > 0)
3434     {
3435         printf("Largest charge group radii for Van der Waals: %5.3f, %5.3f nm\n",
3436                rvdw1,rvdw2);
3437     }
3438     if (rcoul1 > 0)
3439     {
3440         printf("Largest charge group radii for Coulomb:       %5.3f, %5.3f nm\n",
3441                rcoul1,rcoul2);
3442     }
3443
3444     if (ir->rlist > 0)
3445     {
3446         if (rvdw1  + rvdw2  > ir->rlist ||
3447             rcoul1 + rcoul2 > ir->rlist)
3448         {
3449             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);
3450             warning(wi,warn_buf);
3451         }
3452         else
3453         {
3454             /* Here we do not use the zero at cut-off macro,
3455              * since user defined interactions might purposely
3456              * not be zero at the cut-off.
3457              */
3458             if (EVDW_IS_ZERO_AT_CUTOFF(ir->vdwtype) &&
3459                 rvdw1 + rvdw2 > ir->rlist - ir->rvdw)
3460             {
3461                 sprintf(warn_buf,"The sum of the two largest charge group radii (%f) is larger than rlist (%f) - rvdw (%f)\n",
3462                         rvdw1+rvdw2,
3463                         ir->rlist,ir->rvdw);
3464                 if (ir_NVE(ir))
3465                 {
3466                     warning(wi,warn_buf);
3467                 }
3468                 else
3469                 {
3470                     warning_note(wi,warn_buf);
3471                 }
3472             }
3473             if (EEL_IS_ZERO_AT_CUTOFF(ir->coulombtype) &&
3474                 rcoul1 + rcoul2 > ir->rlistlong - ir->rcoulomb)
3475             {
3476                 sprintf(warn_buf,"The sum of the two largest charge group radii (%f) is larger than %s (%f) - rcoulomb (%f)\n",
3477                         rcoul1+rcoul2,
3478                         ir->rlistlong > ir->rlist ? "rlistlong" : "rlist",
3479                         ir->rlistlong,ir->rcoulomb);
3480                 if (ir_NVE(ir))
3481                 {
3482                     warning(wi,warn_buf);
3483                 }
3484                 else
3485                 {
3486                     warning_note(wi,warn_buf);
3487                 }
3488             }
3489         }
3490     }
3491 }