src/gromacs/gmxpreprocess/readir.c

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