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         if (ir->cutoff_scheme != ecutsGROUP)
1154         {
1155             warning_error(wi,"AdresS simulation supports only cutoff-scheme=group");
1156         }
1157         if (!EI_SD(ir->eI))
1158         {
1159             warning_error(wi,"AdresS simulation supports only stochastic dynamics");
1160         }
1161         if (ir->epc != epcNO)
1162         {
1163             warning_error(wi,"AdresS simulation does not support pressure coupling");
1164         }
1165         if (EEL_FULL(ir->coulombtype))
1166         {
1167             warning_error(wi,"AdresS simulation does not support long-range electrostatics");
1168         }
1169     }
1170 }
1171
1172 /* count the number of text elemets separated by whitespace in a string.
1173     str = the input string
1174     maxptr = the maximum number of allowed elements
1175     ptr = the output array of pointers to the first character of each element
1176     returns: the number of elements. */
1177 int str_nelem(const char *str,int maxptr,char *ptr[])
1178 {
1179   int  np=0;
1180   char *copy0,*copy;
1181
1182   copy0=strdup(str);
1183   copy=copy0;
1184   ltrim(copy);
1185   while (*copy != '\0') {
1186     if (np >= maxptr)
1187       gmx_fatal(FARGS,"Too many groups on line: '%s' (max is %d)",
1188                   str,maxptr);
1189     if (ptr)
1190       ptr[np]=copy;
1191     np++;
1192     while ((*copy != '\0') && !isspace(*copy))
1193       copy++;
1194     if (*copy != '\0') {
1195       *copy='\0';
1196       copy++;
1197     }
1198     ltrim(copy);
1199   }
1200   if (ptr == NULL)
1201     sfree(copy0);
1202
1203   return np;
1204 }
1205
1206 /* interpret a number of doubles from a string and put them in an array,
1207    after allocating space for them.
1208    str = the input string
1209    n = the (pre-allocated) number of doubles read
1210    r = the output array of doubles. */
1211 static void parse_n_real(char *str,int *n,real **r)
1212 {
1213   char *ptr[MAXPTR];
1214   int  i;
1215
1216   *n = str_nelem(str,MAXPTR,ptr);
1217
1218   snew(*r,*n);
1219   for(i=0; i<*n; i++) {
1220     (*r)[i] = strtod(ptr[i],NULL);
1221   }
1222 }
1223
1224 static void do_fep_params(t_inputrec *ir, char fep_lambda[][STRLEN],char weights[STRLEN]) {
1225
1226     int i,j,max_n_lambda,nweights,nfep[efptNR];
1227     t_lambda *fep = ir->fepvals;
1228     t_expanded *expand = ir->expandedvals;
1229     real **count_fep_lambdas;
1230     gmx_bool bOneLambda = TRUE;
1231
1232     snew(count_fep_lambdas,efptNR);
1233
1234     /* FEP input processing */
1235     /* first, identify the number of lambda values for each type.
1236        All that are nonzero must have the same number */
1237
1238     for (i=0;i<efptNR;i++)
1239     {
1240         parse_n_real(fep_lambda[i],&(nfep[i]),&(count_fep_lambdas[i]));
1241     }
1242
1243     /* now, determine the number of components.  All must be either zero, or equal. */
1244
1245     max_n_lambda = 0;
1246     for (i=0;i<efptNR;i++)
1247     {
1248         if (nfep[i] > max_n_lambda) {
1249             max_n_lambda = nfep[i];  /* here's a nonzero one.  All of them
1250                                         must have the same number if its not zero.*/
1251             break;
1252         }
1253     }
1254
1255     for (i=0;i<efptNR;i++)
1256     {
1257         if (nfep[i] == 0)
1258         {
1259             ir->fepvals->separate_dvdl[i] = FALSE;
1260         }
1261         else if (nfep[i] == max_n_lambda)
1262         {
1263             if (i!=efptTEMPERATURE)  /* we treat this differently -- not really a reason to compute the derivative with
1264                                         respect to the temperature currently */
1265             {
1266                 ir->fepvals->separate_dvdl[i] = TRUE;
1267             }
1268         }
1269         else
1270         {
1271             gmx_fatal(FARGS,"Number of lambdas (%d) for FEP type %s not equal to number of other types (%d)",
1272                       nfep[i],efpt_names[i],max_n_lambda);
1273         }
1274     }
1275     /* we don't print out dhdl if the temperature is changing, since we can't correctly define dhdl in this case */
1276     ir->fepvals->separate_dvdl[efptTEMPERATURE] = FALSE;
1277
1278     /* the number of lambdas is the number we've read in, which is either zero
1279        or the same for all */
1280     fep->n_lambda = max_n_lambda;
1281
1282     /* allocate space for the array of lambda values */
1283     snew(fep->all_lambda,efptNR);
1284     /* if init_lambda is defined, we need to set lambda */
1285     if ((fep->init_lambda > 0) && (fep->n_lambda == 0))
1286     {
1287         ir->fepvals->separate_dvdl[efptFEP] = TRUE;
1288     }
1289     /* otherwise allocate the space for all of the lambdas, and transfer the data */
1290     for (i=0;i<efptNR;i++)
1291     {
1292         snew(fep->all_lambda[i],fep->n_lambda);
1293         if (nfep[i] > 0)  /* if it's zero, then the count_fep_lambda arrays
1294                              are zero */
1295         {
1296             for (j=0;j<fep->n_lambda;j++)
1297             {
1298                 fep->all_lambda[i][j] = (double)count_fep_lambdas[i][j];
1299             }
1300             sfree(count_fep_lambdas[i]);
1301         }
1302     }
1303     sfree(count_fep_lambdas);
1304
1305     /* "fep-vals" is either zero or the full number. If zero, we'll need to define fep-lambdas for internal
1306        bookkeeping -- for now, init_lambda */
1307
1308     if ((nfep[efptFEP] == 0) && (fep->init_lambda >= 0) && (fep->init_lambda <= 1))
1309     {
1310         for (i=0;i<fep->n_lambda;i++)
1311         {
1312             fep->all_lambda[efptFEP][i] = fep->init_lambda;
1313         }
1314     }
1315
1316     /* check to see if only a single component lambda is defined, and soft core is defined.
1317        In this case, turn on coulomb soft core */
1318
1319     if (max_n_lambda == 0)
1320     {
1321         bOneLambda = TRUE;
1322     }
1323     else
1324     {
1325         for (i=0;i<efptNR;i++)
1326         {
1327             if ((nfep[i] != 0) && (i!=efptFEP))
1328             {
1329                 bOneLambda = FALSE;
1330             }
1331         }
1332     }
1333     if ((bOneLambda) && (fep->sc_alpha > 0))
1334     {
1335         fep->bScCoul = TRUE;
1336     }
1337
1338     /* Fill in the others with the efptFEP if they are not explicitly
1339        specified (i.e. nfep[i] == 0).  This means if fep is not defined,
1340        they are all zero. */
1341
1342     for (i=0;i<efptNR;i++)
1343     {
1344         if ((nfep[i] == 0) && (i!=efptFEP))
1345         {
1346             for (j=0;j<fep->n_lambda;j++)
1347             {
1348                 fep->all_lambda[i][j] = fep->all_lambda[efptFEP][j];
1349             }
1350         }
1351     }
1352
1353
1354     /* make it easier if sc_r_power = 48 by increasing it to the 4th power, to be in the right scale. */
1355     if (fep->sc_r_power == 48)
1356     {
1357         if (fep->sc_alpha > 0.1)
1358         {
1359             gmx_fatal(FARGS,"sc_alpha (%f) for sc_r_power = 48 should usually be between 0.001 and 0.004", fep->sc_alpha);
1360         }
1361     }
1362
1363     expand = ir->expandedvals;
1364     /* now read in the weights */
1365     parse_n_real(weights,&nweights,&(expand->init_lambda_weights));
1366     if (nweights == 0)
1367     {
1368         expand->bInit_weights = FALSE;
1369         snew(expand->init_lambda_weights,fep->n_lambda); /* initialize to zero */
1370     }
1371     else if (nweights != fep->n_lambda)
1372     {
1373         gmx_fatal(FARGS,"Number of weights (%d) is not equal to number of lambda values (%d)",
1374                   nweights,fep->n_lambda);
1375     }
1376     else
1377     {
1378         expand->bInit_weights = TRUE;
1379     }
1380     if ((expand->nstexpanded < 0) && (ir->efep != efepNO)) {
1381         expand->nstexpanded = fep->nstdhdl;
1382         /* if you don't specify nstexpanded when doing expanded ensemble free energy calcs, it is set to nstdhdl */
1383     }
1384     if ((expand->nstexpanded < 0) && ir->bSimTemp) {
1385         expand->nstexpanded = 2*(int)(ir->opts.tau_t[0]/ir->delta_t);
1386         /* if you don't specify nstexpanded when doing expanded ensemble simulated tempering, it is set to
1387            2*tau_t just to be careful so it's not to frequent  */
1388     }
1389 }
1390
1391
1392 static void do_simtemp_params(t_inputrec *ir) {
1393
1394     snew(ir->simtempvals->temperatures,ir->fepvals->n_lambda);
1395     GetSimTemps(ir->fepvals->n_lambda,ir->simtempvals,ir->fepvals->all_lambda[efptTEMPERATURE]);
1396
1397     return;
1398 }
1399
1400 static void do_wall_params(t_inputrec *ir,
1401                            char *wall_atomtype, char *wall_density,
1402                            t_gromppopts *opts)
1403 {
1404     int  nstr,i;
1405     char *names[MAXPTR];
1406     double dbl;
1407
1408     opts->wall_atomtype[0] = NULL;
1409     opts->wall_atomtype[1] = NULL;
1410
1411     ir->wall_atomtype[0] = -1;
1412     ir->wall_atomtype[1] = -1;
1413     ir->wall_density[0] = 0;
1414     ir->wall_density[1] = 0;
1415
1416     if (ir->nwall > 0)
1417     {
1418         nstr = str_nelem(wall_atomtype,MAXPTR,names);
1419         if (nstr != ir->nwall)
1420         {
1421             gmx_fatal(FARGS,"Expected %d elements for wall_atomtype, found %d",
1422                       ir->nwall,nstr);
1423         }
1424         for(i=0; i<ir->nwall; i++)
1425         {
1426             opts->wall_atomtype[i] = strdup(names[i]);
1427         }
1428
1429         if (ir->wall_type == ewt93 || ir->wall_type == ewt104) {
1430             nstr = str_nelem(wall_density,MAXPTR,names);
1431             if (nstr != ir->nwall)
1432             {
1433                 gmx_fatal(FARGS,"Expected %d elements for wall-density, found %d",ir->nwall,nstr);
1434             }
1435             for(i=0; i<ir->nwall; i++)
1436             {
1437                 sscanf(names[i],"%lf",&dbl);
1438                 if (dbl <= 0)
1439                 {
1440                     gmx_fatal(FARGS,"wall-density[%d] = %f\n",i,dbl);
1441                 }
1442                 ir->wall_density[i] = dbl;
1443             }
1444         }
1445     }
1446 }
1447
1448 static void add_wall_energrps(gmx_groups_t *groups,int nwall,t_symtab *symtab)
1449 {
1450   int  i;
1451   t_grps *grps;
1452   char str[STRLEN];
1453
1454   if (nwall > 0) {
1455     srenew(groups->grpname,groups->ngrpname+nwall);
1456     grps = &(groups->grps[egcENER]);
1457     srenew(grps->nm_ind,grps->nr+nwall);
1458     for(i=0; i<nwall; i++) {
1459       sprintf(str,"wall%d",i);
1460       groups->grpname[groups->ngrpname] = put_symtab(symtab,str);
1461       grps->nm_ind[grps->nr++] = groups->ngrpname++;
1462     }
1463   }
1464 }
1465
1466 void read_expandedparams(int *ninp_p,t_inpfile **inp_p,
1467                          t_expanded *expand,warninp_t wi)
1468 {
1469   int  ninp,nerror=0;
1470   t_inpfile *inp;
1471
1472   ninp   = *ninp_p;
1473   inp    = *inp_p;
1474
1475   /* read expanded ensemble parameters */
1476   CCTYPE ("expanded ensemble variables");
1477   ITYPE ("nstexpanded",expand->nstexpanded,-1);
1478   EETYPE("lmc-stats", expand->elamstats, elamstats_names);
1479   EETYPE("lmc-move", expand->elmcmove, elmcmove_names);
1480   EETYPE("lmc-weights-equil",expand->elmceq,elmceq_names);
1481   ITYPE ("weight-equil-number-all-lambda",expand->equil_n_at_lam,-1);
1482   ITYPE ("weight-equil-number-samples",expand->equil_samples,-1);
1483   ITYPE ("weight-equil-number-steps",expand->equil_steps,-1);
1484   RTYPE ("weight-equil-wl-delta",expand->equil_wl_delta,-1);
1485   RTYPE ("weight-equil-count-ratio",expand->equil_ratio,-1);
1486   CCTYPE("Seed for Monte Carlo in lambda space");
1487   ITYPE ("lmc-seed",expand->lmc_seed,-1);
1488   RTYPE ("mc-temperature",expand->mc_temp,-1);
1489   ITYPE ("lmc-repeats",expand->lmc_repeats,1);
1490   ITYPE ("lmc-gibbsdelta",expand->gibbsdeltalam,-1);
1491   ITYPE ("lmc-forced-nstart",expand->lmc_forced_nstart,0);
1492   EETYPE("symmetrized-transition-matrix", expand->bSymmetrizedTMatrix, yesno_names);
1493   ITYPE("nst-transition-matrix", expand->nstTij, -1);
1494   ITYPE ("mininum-var-min",expand->minvarmin, 100); /*default is reasonable */
1495   ITYPE ("weight-c-range",expand->c_range, 0); /* default is just C=0 */
1496   RTYPE ("wl-scale",expand->wl_scale,0.8);
1497   RTYPE ("wl-ratio",expand->wl_ratio,0.8);
1498   RTYPE ("init-wl-delta",expand->init_wl_delta,1.0);
1499   EETYPE("wl-oneovert",expand->bWLoneovert,yesno_names);
1500
1501   *ninp_p   = ninp;
1502   *inp_p    = inp;
1503
1504   return;
1505 }
1506
1507 void get_ir(const char *mdparin,const char *mdparout,
1508             t_inputrec *ir,t_gromppopts *opts,
1509             warninp_t wi)
1510 {
1511   char      *dumstr[2];
1512   double    dumdub[2][6];
1513   t_inpfile *inp;
1514   const char *tmp;
1515   int       i,j,m,ninp;
1516   char      warn_buf[STRLEN];
1517   t_lambda  *fep = ir->fepvals;
1518   t_expanded *expand = ir->expandedvals;
1519
1520   inp = read_inpfile(mdparin, &ninp, NULL, wi);
1521
1522   snew(dumstr[0],STRLEN);
1523   snew(dumstr[1],STRLEN);
1524
1525   /* remove the following deprecated commands */
1526   REM_TYPE("title");
1527   REM_TYPE("cpp");
1528   REM_TYPE("domain-decomposition");
1529   REM_TYPE("andersen-seed");
1530   REM_TYPE("dihre");
1531   REM_TYPE("dihre-fc");
1532   REM_TYPE("dihre-tau");
1533   REM_TYPE("nstdihreout");
1534   REM_TYPE("nstcheckpoint");
1535
1536   /* replace the following commands with the clearer new versions*/
1537   REPL_TYPE("unconstrained-start","continuation");
1538   REPL_TYPE("foreign-lambda","fep-lambdas");
1539
1540   CCTYPE ("VARIOUS PREPROCESSING OPTIONS");
1541   CTYPE ("Preprocessor information: use cpp syntax.");
1542   CTYPE ("e.g.: -I/home/joe/doe -I/home/mary/roe");
1543   STYPE ("include",     opts->include,  NULL);
1544   CTYPE ("e.g.: -DPOSRES -DFLEXIBLE (note these variable names are case sensitive)");
1545   STYPE ("define",      opts->define,   NULL);
1546
1547   CCTYPE ("RUN CONTROL PARAMETERS");
1548   EETYPE("integrator",  ir->eI,         ei_names);
1549   CTYPE ("Start time and timestep in ps");
1550   RTYPE ("tinit",       ir->init_t,     0.0);
1551   RTYPE ("dt",          ir->delta_t,    0.001);
1552   STEPTYPE ("nsteps",   ir->nsteps,     0);
1553   CTYPE ("For exact run continuation or redoing part of a run");
1554   STEPTYPE ("init-step",ir->init_step,  0);
1555   CTYPE ("Part index is updated automatically on checkpointing (keeps files separate)");
1556   ITYPE ("simulation-part", ir->simulation_part, 1);
1557   CTYPE ("mode for center of mass motion removal");
1558   EETYPE("comm-mode",   ir->comm_mode,  ecm_names);
1559   CTYPE ("number of steps for center of mass motion removal");
1560   ITYPE ("nstcomm",     ir->nstcomm,    100);
1561   CTYPE ("group(s) for center of mass motion removal");
1562   STYPE ("comm-grps",   vcm,            NULL);
1563
1564   CCTYPE ("LANGEVIN DYNAMICS OPTIONS");
1565   CTYPE ("Friction coefficient (amu/ps) and random seed");
1566   RTYPE ("bd-fric",     ir->bd_fric,    0.0);
1567   ITYPE ("ld-seed",     ir->ld_seed,    1993);
1568
1569   /* Em stuff */
1570   CCTYPE ("ENERGY MINIMIZATION OPTIONS");
1571   CTYPE ("Force tolerance and initial step-size");
1572   RTYPE ("emtol",       ir->em_tol,     10.0);
1573   RTYPE ("emstep",      ir->em_stepsize,0.01);
1574   CTYPE ("Max number of iterations in relax-shells");
1575   ITYPE ("niter",       ir->niter,      20);
1576   CTYPE ("Step size (ps^2) for minimization of flexible constraints");
1577   RTYPE ("fcstep",      ir->fc_stepsize, 0);
1578   CTYPE ("Frequency of steepest descents steps when doing CG");
1579   ITYPE ("nstcgsteep",  ir->nstcgsteep, 1000);
1580   ITYPE ("nbfgscorr",   ir->nbfgscorr,  10);
1581
1582   CCTYPE ("TEST PARTICLE INSERTION OPTIONS");
1583   RTYPE ("rtpi",        ir->rtpi,       0.05);
1584
1585   /* Output options */
1586   CCTYPE ("OUTPUT CONTROL OPTIONS");
1587   CTYPE ("Output frequency for coords (x), velocities (v) and forces (f)");
1588   ITYPE ("nstxout",     ir->nstxout,    0);
1589   ITYPE ("nstvout",     ir->nstvout,    0);
1590   ITYPE ("nstfout",     ir->nstfout,    0);
1591   ir->nstcheckpoint = 1000;
1592   CTYPE ("Output frequency for energies to log file and energy file");
1593   ITYPE ("nstlog",      ir->nstlog,     1000);
1594   ITYPE ("nstcalcenergy",ir->nstcalcenergy,     100);
1595   ITYPE ("nstenergy",   ir->nstenergy,  1000);
1596   CTYPE ("Output frequency and precision for .xtc file");
1597   ITYPE ("nstxtcout",   ir->nstxtcout,  0);
1598   RTYPE ("xtc-precision",ir->xtcprec,   1000.0);
1599   CTYPE ("This selects the subset of atoms for the .xtc file. You can");
1600   CTYPE ("select multiple groups. By default all atoms will be written.");
1601   STYPE ("xtc-grps",    xtc_grps,       NULL);
1602   CTYPE ("Selection of energy groups");
1603   STYPE ("energygrps",  energy,         NULL);
1604
1605   /* Neighbor searching */
1606   CCTYPE ("NEIGHBORSEARCHING PARAMETERS");
1607   CTYPE ("cut-off scheme (group: using charge groups, Verlet: particle based cut-offs)");
1608   EETYPE("cutoff-scheme",     ir->cutoff_scheme,    ecutscheme_names);
1609   CTYPE ("nblist update frequency");
1610   ITYPE ("nstlist",     ir->nstlist,    10);
1611   CTYPE ("ns algorithm (simple or grid)");
1612   EETYPE("ns-type",     ir->ns_type,    ens_names);
1613   /* set ndelta to the optimal value of 2 */
1614   ir->ndelta = 2;
1615   CTYPE ("Periodic boundary conditions: xyz, no, xy");
1616   EETYPE("pbc",         ir->ePBC,       epbc_names);
1617   EETYPE("periodic-molecules", ir->bPeriodicMols, yesno_names);
1618   CTYPE ("Allowed energy drift due to the Verlet buffer in kJ/mol/ps per atom,");
1619   CTYPE ("a value of -1 means: use rlist");
1620   RTYPE("verlet-buffer-drift", ir->verletbuf_drift,    0.005);
1621   CTYPE ("nblist cut-off");
1622   RTYPE ("rlist",       ir->rlist,      1.0);
1623   CTYPE ("long-range cut-off for switched potentials");
1624   RTYPE ("rlistlong",   ir->rlistlong,  -1);
1625   ITYPE ("nstcalclr",   ir->nstcalclr,  -1);
1626
1627   /* Electrostatics */
1628   CCTYPE ("OPTIONS FOR ELECTROSTATICS AND VDW");
1629   CTYPE ("Method for doing electrostatics");
1630   EETYPE("coulombtype", ir->coulombtype,    eel_names);
1631   EETYPE("coulomb-modifier",    ir->coulomb_modifier,    eintmod_names);
1632   CTYPE ("cut-off lengths");
1633   RTYPE ("rcoulomb-switch",     ir->rcoulomb_switch,    0.0);
1634   RTYPE ("rcoulomb",    ir->rcoulomb,   1.0);
1635   CTYPE ("Relative dielectric constant for the medium and the reaction field");
1636   RTYPE ("epsilon-r",   ir->epsilon_r,  1.0);
1637   RTYPE ("epsilon-rf",  ir->epsilon_rf, 0.0);
1638   CTYPE ("Method for doing Van der Waals");
1639   EETYPE("vdw-type",    ir->vdwtype,    evdw_names);
1640   EETYPE("vdw-modifier",        ir->vdw_modifier,    eintmod_names);
1641   CTYPE ("cut-off lengths");
1642   RTYPE ("rvdw-switch", ir->rvdw_switch,        0.0);
1643   RTYPE ("rvdw",        ir->rvdw,       1.0);
1644   CTYPE ("Apply long range dispersion corrections for Energy and Pressure");
1645   EETYPE("DispCorr",    ir->eDispCorr,  edispc_names);
1646   CTYPE ("Extension of the potential lookup tables beyond the cut-off");
1647   RTYPE ("table-extension", ir->tabext, 1.0);
1648   CTYPE ("Separate tables between energy group pairs");
1649   STYPE ("energygrp-table", egptable,   NULL);
1650   CTYPE ("Spacing for the PME/PPPM FFT grid");
1651   RTYPE ("fourierspacing", ir->fourier_spacing,0.12);
1652   CTYPE ("FFT grid size, when a value is 0 fourierspacing will be used");
1653   ITYPE ("fourier-nx",  ir->nkx,         0);
1654   ITYPE ("fourier-ny",  ir->nky,         0);
1655   ITYPE ("fourier-nz",  ir->nkz,         0);
1656   CTYPE ("EWALD/PME/PPPM parameters");
1657   ITYPE ("pme-order",   ir->pme_order,   4);
1658   RTYPE ("ewald-rtol",  ir->ewald_rtol, 0.00001);
1659   EETYPE("ewald-geometry", ir->ewald_geometry, eewg_names);
1660   RTYPE ("epsilon-surface", ir->epsilon_surface, 0.0);
1661   EETYPE("optimize-fft",ir->bOptFFT,  yesno_names);
1662
1663   CCTYPE("IMPLICIT SOLVENT ALGORITHM");
1664   EETYPE("implicit-solvent", ir->implicit_solvent, eis_names);
1665
1666   CCTYPE ("GENERALIZED BORN ELECTROSTATICS");
1667   CTYPE ("Algorithm for calculating Born radii");
1668   EETYPE("gb-algorithm", ir->gb_algorithm, egb_names);
1669   CTYPE ("Frequency of calculating the Born radii inside rlist");
1670   ITYPE ("nstgbradii", ir->nstgbradii, 1);
1671   CTYPE ("Cutoff for Born radii calculation; the contribution from atoms");
1672   CTYPE ("between rlist and rgbradii is updated every nstlist steps");
1673   RTYPE ("rgbradii",  ir->rgbradii, 1.0);
1674   CTYPE ("Dielectric coefficient of the implicit solvent");
1675   RTYPE ("gb-epsilon-solvent",ir->gb_epsilon_solvent, 80.0);
1676   CTYPE ("Salt concentration in M for Generalized Born models");
1677   RTYPE ("gb-saltconc",  ir->gb_saltconc, 0.0);
1678   CTYPE ("Scaling factors used in the OBC GB model. Default values are OBC(II)");
1679   RTYPE ("gb-obc-alpha", ir->gb_obc_alpha, 1.0);
1680   RTYPE ("gb-obc-beta", ir->gb_obc_beta, 0.8);
1681   RTYPE ("gb-obc-gamma", ir->gb_obc_gamma, 4.85);
1682   RTYPE ("gb-dielectric-offset", ir->gb_dielectric_offset, 0.009);
1683   EETYPE("sa-algorithm", ir->sa_algorithm, esa_names);
1684   CTYPE ("Surface tension (kJ/mol/nm^2) for the SA (nonpolar surface) part of GBSA");
1685   CTYPE ("The value -1 will set default value for Still/HCT/OBC GB-models.");
1686   RTYPE ("sa-surface-tension", ir->sa_surface_tension, -1);
1687
1688   /* Coupling stuff */
1689   CCTYPE ("OPTIONS FOR WEAK COUPLING ALGORITHMS");
1690   CTYPE ("Temperature coupling");
1691   EETYPE("tcoupl",      ir->etc,        etcoupl_names);
1692   ITYPE ("nsttcouple", ir->nsttcouple,  -1);
1693   ITYPE("nh-chain-length",     ir->opts.nhchainlength, NHCHAINLENGTH);
1694   EETYPE("print-nose-hoover-chain-variables", ir->bPrintNHChains, yesno_names);
1695   CTYPE ("Groups to couple separately");
1696   STYPE ("tc-grps",     tcgrps,         NULL);
1697   CTYPE ("Time constant (ps) and reference temperature (K)");
1698   STYPE ("tau-t",       tau_t,          NULL);
1699   STYPE ("ref-t",       ref_t,          NULL);
1700   CTYPE ("pressure coupling");
1701   EETYPE("pcoupl",      ir->epc,        epcoupl_names);
1702   EETYPE("pcoupltype",  ir->epct,       epcoupltype_names);
1703   ITYPE ("nstpcouple", ir->nstpcouple,  -1);
1704   CTYPE ("Time constant (ps), compressibility (1/bar) and reference P (bar)");
1705   RTYPE ("tau-p",       ir->tau_p,      1.0);
1706   STYPE ("compressibility",     dumstr[0],      NULL);
1707   STYPE ("ref-p",       dumstr[1],      NULL);
1708   CTYPE ("Scaling of reference coordinates, No, All or COM");
1709   EETYPE ("refcoord-scaling",ir->refcoord_scaling,erefscaling_names);
1710
1711   /* QMMM */
1712   CCTYPE ("OPTIONS FOR QMMM calculations");
1713   EETYPE("QMMM", ir->bQMMM, yesno_names);
1714   CTYPE ("Groups treated Quantum Mechanically");
1715   STYPE ("QMMM-grps",  QMMM,          NULL);
1716   CTYPE ("QM method");
1717   STYPE("QMmethod",     QMmethod, NULL);
1718   CTYPE ("QMMM scheme");
1719   EETYPE("QMMMscheme",  ir->QMMMscheme,    eQMMMscheme_names);
1720   CTYPE ("QM basisset");
1721   STYPE("QMbasis",      QMbasis, NULL);
1722   CTYPE ("QM charge");
1723   STYPE ("QMcharge",    QMcharge,NULL);
1724   CTYPE ("QM multiplicity");
1725   STYPE ("QMmult",      QMmult,NULL);
1726   CTYPE ("Surface Hopping");
1727   STYPE ("SH",          bSH, NULL);
1728   CTYPE ("CAS space options");
1729   STYPE ("CASorbitals",      CASorbitals,   NULL);
1730   STYPE ("CASelectrons",     CASelectrons,  NULL);
1731   STYPE ("SAon", SAon, NULL);
1732   STYPE ("SAoff",SAoff,NULL);
1733   STYPE ("SAsteps",  SAsteps, NULL);
1734   CTYPE ("Scale factor for MM charges");
1735   RTYPE ("MMChargeScaleFactor", ir->scalefactor, 1.0);
1736   CTYPE ("Optimization of QM subsystem");
1737   STYPE ("bOPT",          bOPT, NULL);
1738   STYPE ("bTS",          bTS, NULL);
1739
1740   /* Simulated annealing */
1741   CCTYPE("SIMULATED ANNEALING");
1742   CTYPE ("Type of annealing for each temperature group (no/single/periodic)");
1743   STYPE ("annealing",   anneal,      NULL);
1744   CTYPE ("Number of time points to use for specifying annealing in each group");
1745   STYPE ("annealing-npoints", anneal_npoints, NULL);
1746   CTYPE ("List of times at the annealing points for each group");
1747   STYPE ("annealing-time",       anneal_time,       NULL);
1748   CTYPE ("Temp. at each annealing point, for each group.");
1749   STYPE ("annealing-temp",  anneal_temp,  NULL);
1750
1751   /* Startup run */
1752   CCTYPE ("GENERATE VELOCITIES FOR STARTUP RUN");
1753   EETYPE("gen-vel",     opts->bGenVel,  yesno_names);
1754   RTYPE ("gen-temp",    opts->tempi,    300.0);
1755   ITYPE ("gen-seed",    opts->seed,     173529);
1756
1757   /* Shake stuff */
1758   CCTYPE ("OPTIONS FOR BONDS");
1759   EETYPE("constraints", opts->nshake,   constraints);
1760   CTYPE ("Type of constraint algorithm");
1761   EETYPE("constraint-algorithm",  ir->eConstrAlg, econstr_names);
1762   CTYPE ("Do not constrain the start configuration");
1763   EETYPE("continuation", ir->bContinuation, yesno_names);
1764   CTYPE ("Use successive overrelaxation to reduce the number of shake iterations");
1765   EETYPE("Shake-SOR", ir->bShakeSOR, yesno_names);
1766   CTYPE ("Relative tolerance of shake");
1767   RTYPE ("shake-tol", ir->shake_tol, 0.0001);
1768   CTYPE ("Highest order in the expansion of the constraint coupling matrix");
1769   ITYPE ("lincs-order", ir->nProjOrder, 4);
1770   CTYPE ("Number of iterations in the final step of LINCS. 1 is fine for");
1771   CTYPE ("normal simulations, but use 2 to conserve energy in NVE runs.");
1772   CTYPE ("For energy minimization with constraints it should be 4 to 8.");
1773   ITYPE ("lincs-iter", ir->nLincsIter, 1);
1774   CTYPE ("Lincs will write a warning to the stderr if in one step a bond");
1775   CTYPE ("rotates over more degrees than");
1776   RTYPE ("lincs-warnangle", ir->LincsWarnAngle, 30.0);
1777   CTYPE ("Convert harmonic bonds to morse potentials");
1778   EETYPE("morse",       opts->bMorse,yesno_names);
1779
1780   /* Energy group exclusions */
1781   CCTYPE ("ENERGY GROUP EXCLUSIONS");
1782   CTYPE ("Pairs of energy groups for which all non-bonded interactions are excluded");
1783   STYPE ("energygrp-excl", egpexcl,     NULL);
1784
1785   /* Walls */
1786   CCTYPE ("WALLS");
1787   CTYPE ("Number of walls, type, atom types, densities and box-z scale factor for Ewald");
1788   ITYPE ("nwall", ir->nwall, 0);
1789   EETYPE("wall-type",     ir->wall_type,   ewt_names);
1790   RTYPE ("wall-r-linpot", ir->wall_r_linpot, -1);
1791   STYPE ("wall-atomtype", wall_atomtype, NULL);
1792   STYPE ("wall-density",  wall_density,  NULL);
1793   RTYPE ("wall-ewald-zfac", ir->wall_ewald_zfac, 3);
1794
1795   /* COM pulling */
1796   CCTYPE("COM PULLING");
1797   CTYPE("Pull type: no, umbrella, constraint or constant-force");
1798   EETYPE("pull",          ir->ePull, epull_names);
1799   if (ir->ePull != epullNO) {
1800     snew(ir->pull,1);
1801     pull_grp = read_pullparams(&ninp,&inp,ir->pull,&opts->pull_start,wi);
1802   }
1803
1804   /* Enforced rotation */
1805   CCTYPE("ENFORCED ROTATION");
1806   CTYPE("Enforced rotation: No or Yes");
1807   EETYPE("rotation",       ir->bRot, yesno_names);
1808   if (ir->bRot) {
1809     snew(ir->rot,1);
1810     rot_grp = read_rotparams(&ninp,&inp,ir->rot,wi);
1811   }
1812
1813   /* Refinement */
1814   CCTYPE("NMR refinement stuff");
1815   CTYPE ("Distance restraints type: No, Simple or Ensemble");
1816   EETYPE("disre",       ir->eDisre,     edisre_names);
1817   CTYPE ("Force weighting of pairs in one distance restraint: Conservative or Equal");
1818   EETYPE("disre-weighting", ir->eDisreWeighting, edisreweighting_names);
1819   CTYPE ("Use sqrt of the time averaged times the instantaneous violation");
1820   EETYPE("disre-mixed", ir->bDisreMixed, yesno_names);
1821   RTYPE ("disre-fc",    ir->dr_fc,      1000.0);
1822   RTYPE ("disre-tau",   ir->dr_tau,     0.0);
1823   CTYPE ("Output frequency for pair distances to energy file");
1824   ITYPE ("nstdisreout", ir->nstdisreout, 100);
1825   CTYPE ("Orientation restraints: No or Yes");
1826   EETYPE("orire",       opts->bOrire,   yesno_names);
1827   CTYPE ("Orientation restraints force constant and tau for time averaging");
1828   RTYPE ("orire-fc",    ir->orires_fc,  0.0);
1829   RTYPE ("orire-tau",   ir->orires_tau, 0.0);
1830   STYPE ("orire-fitgrp",orirefitgrp,    NULL);
1831   CTYPE ("Output frequency for trace(SD) and S to energy file");
1832   ITYPE ("nstorireout", ir->nstorireout, 100);
1833
1834   /* free energy variables */
1835   CCTYPE ("Free energy variables");
1836   EETYPE("free-energy", ir->efep, efep_names);
1837   STYPE ("couple-moltype",  couple_moltype,  NULL);
1838   EETYPE("couple-lambda0", opts->couple_lam0, couple_lam);
1839   EETYPE("couple-lambda1", opts->couple_lam1, couple_lam);
1840   EETYPE("couple-intramol", opts->bCoupleIntra, yesno_names);
1841
1842   RTYPE ("init-lambda", fep->init_lambda,-1); /* start with -1 so
1843                                                  we can recognize if
1844                                                  it was not entered */
1845   ITYPE ("init-lambda-state", fep->init_fep_state,0);
1846   RTYPE ("delta-lambda",fep->delta_lambda,0.0);
1847   ITYPE ("nstdhdl",fep->nstdhdl, 100);
1848   STYPE ("fep-lambdas", fep_lambda[efptFEP], NULL);
1849   STYPE ("mass-lambdas", fep_lambda[efptMASS], NULL);
1850   STYPE ("coul-lambdas", fep_lambda[efptCOUL], NULL);
1851   STYPE ("vdw-lambdas", fep_lambda[efptVDW], NULL);
1852   STYPE ("bonded-lambdas", fep_lambda[efptBONDED], NULL);
1853   STYPE ("restraint-lambdas", fep_lambda[efptRESTRAINT], NULL);
1854   STYPE ("temperature-lambdas", fep_lambda[efptTEMPERATURE], NULL);
1855   STYPE ("init-lambda-weights",lambda_weights,NULL);
1856   EETYPE("dhdl-print-energy", fep->bPrintEnergy, yesno_names);
1857   RTYPE ("sc-alpha",fep->sc_alpha,0.0);
1858   ITYPE ("sc-power",fep->sc_power,1);
1859   RTYPE ("sc-r-power",fep->sc_r_power,6.0);
1860   RTYPE ("sc-sigma",fep->sc_sigma,0.3);
1861   EETYPE("sc-coul",fep->bScCoul,yesno_names);
1862   ITYPE ("dh_hist_size", fep->dh_hist_size, 0);
1863   RTYPE ("dh_hist_spacing", fep->dh_hist_spacing, 0.1);
1864   EETYPE("separate-dhdl-file", fep->separate_dhdl_file,
1865                                separate_dhdl_file_names);
1866   EETYPE("dhdl-derivatives", fep->dhdl_derivatives, dhdl_derivatives_names);
1867   ITYPE ("dh_hist_size", fep->dh_hist_size, 0);
1868   RTYPE ("dh_hist_spacing", fep->dh_hist_spacing, 0.1);
1869
1870   /* Non-equilibrium MD stuff */
1871   CCTYPE("Non-equilibrium MD stuff");
1872   STYPE ("acc-grps",    accgrps,        NULL);
1873   STYPE ("accelerate",  acc,            NULL);
1874   STYPE ("freezegrps",  freeze,         NULL);
1875   STYPE ("freezedim",   frdim,          NULL);
1876   RTYPE ("cos-acceleration", ir->cos_accel, 0);
1877   STYPE ("deform",      deform,         NULL);
1878
1879   /* simulated tempering variables */
1880   CCTYPE("simulated tempering variables");
1881   EETYPE("simulated-tempering",ir->bSimTemp,yesno_names);
1882   EETYPE("simulated-tempering-scaling",ir->simtempvals->eSimTempScale,esimtemp_names);
1883   RTYPE("sim-temp-low",ir->simtempvals->simtemp_low,300.0);
1884   RTYPE("sim-temp-high",ir->simtempvals->simtemp_high,300.0);
1885
1886   /* expanded ensemble variables */
1887   if (ir->efep==efepEXPANDED || ir->bSimTemp)
1888   {
1889       read_expandedparams(&ninp,&inp,expand,wi);
1890   }
1891
1892   /* Electric fields */
1893   CCTYPE("Electric fields");
1894   CTYPE ("Format is number of terms (int) and for all terms an amplitude (real)");
1895   CTYPE ("and a phase angle (real)");
1896   STYPE ("E-x",         efield_x,       NULL);
1897   STYPE ("E-xt",        efield_xt,      NULL);
1898   STYPE ("E-y",         efield_y,       NULL);
1899   STYPE ("E-yt",        efield_yt,      NULL);
1900   STYPE ("E-z",         efield_z,       NULL);
1901   STYPE ("E-zt",        efield_zt,      NULL);
1902
1903   /* AdResS defined thingies */
1904   CCTYPE ("AdResS parameters");
1905   EETYPE("adress",       ir->bAdress, yesno_names);
1906   if (ir->bAdress) {
1907     snew(ir->adress,1);
1908     read_adressparams(&ninp,&inp,ir->adress,wi);
1909   }
1910
1911   /* User defined thingies */
1912   CCTYPE ("User defined thingies");
1913   STYPE ("user1-grps",  user1,          NULL);
1914   STYPE ("user2-grps",  user2,          NULL);
1915   ITYPE ("userint1",    ir->userint1,   0);
1916   ITYPE ("userint2",    ir->userint2,   0);
1917   ITYPE ("userint3",    ir->userint3,   0);
1918   ITYPE ("userint4",    ir->userint4,   0);
1919   RTYPE ("userreal1",   ir->userreal1,  0);
1920   RTYPE ("userreal2",   ir->userreal2,  0);
1921   RTYPE ("userreal3",   ir->userreal3,  0);
1922   RTYPE ("userreal4",   ir->userreal4,  0);
1923 #undef CTYPE
1924
1925   write_inpfile(mdparout,ninp,inp,FALSE,wi);
1926   for (i=0; (i<ninp); i++) {
1927     sfree(inp[i].name);
1928     sfree(inp[i].value);
1929   }
1930   sfree(inp);
1931
1932   /* Process options if necessary */
1933   for(m=0; m<2; m++) {
1934     for(i=0; i<2*DIM; i++)
1935       dumdub[m][i]=0.0;
1936     if(ir->epc) {
1937       switch (ir->epct) {
1938       case epctISOTROPIC:
1939         if (sscanf(dumstr[m],"%lf",&(dumdub[m][XX]))!=1) {
1940         warning_error(wi,"Pressure coupling not enough values (I need 1)");
1941         }
1942         dumdub[m][YY]=dumdub[m][ZZ]=dumdub[m][XX];
1943         break;
1944       case epctSEMIISOTROPIC:
1945       case epctSURFACETENSION:
1946         if (sscanf(dumstr[m],"%lf%lf",
1947                    &(dumdub[m][XX]),&(dumdub[m][ZZ]))!=2) {
1948         warning_error(wi,"Pressure coupling not enough values (I need 2)");
1949         }
1950         dumdub[m][YY]=dumdub[m][XX];
1951         break;
1952       case epctANISOTROPIC:
1953         if (sscanf(dumstr[m],"%lf%lf%lf%lf%lf%lf",
1954                    &(dumdub[m][XX]),&(dumdub[m][YY]),&(dumdub[m][ZZ]),
1955                    &(dumdub[m][3]),&(dumdub[m][4]),&(dumdub[m][5]))!=6) {
1956         warning_error(wi,"Pressure coupling not enough values (I need 6)");
1957         }
1958         break;
1959       default:
1960         gmx_fatal(FARGS,"Pressure coupling type %s not implemented yet",
1961                     epcoupltype_names[ir->epct]);
1962       }
1963     }
1964   }
1965   clear_mat(ir->ref_p);
1966   clear_mat(ir->compress);
1967   for(i=0; i<DIM; i++) {
1968     ir->ref_p[i][i]    = dumdub[1][i];
1969     ir->compress[i][i] = dumdub[0][i];
1970   }
1971   if (ir->epct == epctANISOTROPIC) {
1972     ir->ref_p[XX][YY] = dumdub[1][3];
1973     ir->ref_p[XX][ZZ] = dumdub[1][4];
1974     ir->ref_p[YY][ZZ] = dumdub[1][5];
1975     if (ir->ref_p[XX][YY]!=0 && ir->ref_p[XX][ZZ]!=0 && ir->ref_p[YY][ZZ]!=0) {
1976       warning(wi,"All off-diagonal reference pressures are non-zero. Are you sure you want to apply a threefold shear stress?\n");
1977     }
1978     ir->compress[XX][YY] = dumdub[0][3];
1979     ir->compress[XX][ZZ] = dumdub[0][4];
1980     ir->compress[YY][ZZ] = dumdub[0][5];
1981     for(i=0; i<DIM; i++) {
1982       for(m=0; m<i; m++) {
1983         ir->ref_p[i][m] = ir->ref_p[m][i];
1984         ir->compress[i][m] = ir->compress[m][i];
1985       }
1986     }
1987   }
1988
1989   if (ir->comm_mode == ecmNO)
1990     ir->nstcomm = 0;
1991
1992   opts->couple_moltype = NULL;
1993   if (strlen(couple_moltype) > 0)
1994   {
1995       if (ir->efep != efepNO)
1996       {
1997           opts->couple_moltype = strdup(couple_moltype);
1998           if (opts->couple_lam0 == opts->couple_lam1)
1999           {
2000               warning(wi,"The lambda=0 and lambda=1 states for coupling are identical");
2001           }
2002           if (ir->eI == eiMD && (opts->couple_lam0 == ecouplamNONE ||
2003                                  opts->couple_lam1 == ecouplamNONE))
2004           {
2005               warning(wi,"For proper sampling of the (nearly) decoupled state, stochastic dynamics should be used");
2006           }
2007       }
2008       else
2009       {
2010           warning(wi,"Can not couple a molecule with free_energy = no");
2011       }
2012   }
2013   /* FREE ENERGY AND EXPANDED ENSEMBLE OPTIONS */
2014   if (ir->efep != efepNO) {
2015       if (fep->delta_lambda > 0) {
2016           ir->efep = efepSLOWGROWTH;
2017       }
2018   }
2019
2020   if (ir->bSimTemp) {
2021       fep->bPrintEnergy = TRUE;
2022       /* always print out the energy to dhdl if we are doing expanded ensemble, since we need the total energy
2023          if the temperature is changing. */
2024   }
2025
2026   if ((ir->efep != efepNO) || ir->bSimTemp)
2027   {
2028       ir->bExpanded = FALSE;
2029       if ((ir->efep == efepEXPANDED) || ir->bSimTemp)
2030       {
2031           ir->bExpanded = TRUE;
2032       }
2033       do_fep_params(ir,fep_lambda,lambda_weights);
2034       if (ir->bSimTemp) { /* done after fep params */
2035           do_simtemp_params(ir);
2036       }
2037   }
2038   else
2039   {
2040       ir->fepvals->n_lambda = 0;
2041   }
2042
2043   /* WALL PARAMETERS */
2044
2045   do_wall_params(ir,wall_atomtype,wall_density,opts);
2046
2047   /* ORIENTATION RESTRAINT PARAMETERS */
2048
2049   if (opts->bOrire && str_nelem(orirefitgrp,MAXPTR,NULL)!=1) {
2050       warning_error(wi,"ERROR: Need one orientation restraint fit group\n");
2051   }
2052
2053   /* DEFORMATION PARAMETERS */
2054
2055   clear_mat(ir->deform);
2056   for(i=0; i<6; i++)
2057   {
2058       dumdub[0][i] = 0;
2059   }
2060   m = sscanf(deform,"%lf %lf %lf %lf %lf %lf",
2061              &(dumdub[0][0]),&(dumdub[0][1]),&(dumdub[0][2]),
2062              &(dumdub[0][3]),&(dumdub[0][4]),&(dumdub[0][5]));
2063   for(i=0; i<3; i++)
2064   {
2065       ir->deform[i][i] = dumdub[0][i];
2066   }
2067   ir->deform[YY][XX] = dumdub[0][3];
2068   ir->deform[ZZ][XX] = dumdub[0][4];
2069   ir->deform[ZZ][YY] = dumdub[0][5];
2070   if (ir->epc != epcNO) {
2071     for(i=0; i<3; i++)
2072       for(j=0; j<=i; j++)
2073         if (ir->deform[i][j]!=0 && ir->compress[i][j]!=0) {
2074         warning_error(wi,"A box element has deform set and compressibility > 0");
2075         }
2076     for(i=0; i<3; i++)
2077       for(j=0; j<i; j++)
2078         if (ir->deform[i][j]!=0) {
2079           for(m=j; m<DIM; m++)
2080             if (ir->compress[m][j]!=0) {
2081               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.");
2082               warning(wi,warn_buf);
2083             }
2084         }
2085   }
2086
2087   sfree(dumstr[0]);
2088   sfree(dumstr[1]);
2089 }
2090
2091 static int search_QMstring(char *s,int ng,const char *gn[])
2092 {
2093   /* same as normal search_string, but this one searches QM strings */
2094   int i;
2095
2096   for(i=0; (i<ng); i++)
2097     if (gmx_strcasecmp(s,gn[i]) == 0)
2098       return i;
2099
2100   gmx_fatal(FARGS,"this QM method or basisset (%s) is not implemented\n!",s);
2101
2102   return -1;
2103
2104 } /* search_QMstring */
2105
2106
2107 int search_string(char *s,int ng,char *gn[])
2108 {
2109   int i;
2110
2111   for(i=0; (i<ng); i++)
2112   {
2113     if (gmx_strcasecmp(s,gn[i]) == 0)
2114     {
2115       return i;
2116     }
2117   }
2118
2119   gmx_fatal(FARGS,
2120             "Group %s referenced in the .mdp file was not found in the index file.\n"
2121             "Group names must match either [moleculetype] names or custom index group\n"
2122             "names, in which case you must supply an index file to the '-n' option\n"
2123             "of grompp.",
2124             s);
2125
2126   return -1;
2127 }
2128
2129 static gmx_bool do_numbering(int natoms,gmx_groups_t *groups,int ng,char *ptrs[],
2130                          t_blocka *block,char *gnames[],
2131                          int gtype,int restnm,
2132                          int grptp,gmx_bool bVerbose,
2133                          warninp_t wi)
2134 {
2135     unsigned short *cbuf;
2136     t_grps *grps=&(groups->grps[gtype]);
2137     int    i,j,gid,aj,ognr,ntot=0;
2138     const char *title;
2139     gmx_bool   bRest;
2140     char   warn_buf[STRLEN];
2141
2142     if (debug)
2143     {
2144         fprintf(debug,"Starting numbering %d groups of type %d\n",ng,gtype);
2145     }
2146
2147     title = gtypes[gtype];
2148
2149     snew(cbuf,natoms);
2150     /* Mark all id's as not set */
2151     for(i=0; (i<natoms); i++)
2152     {
2153         cbuf[i] = NOGID;
2154     }
2155
2156     snew(grps->nm_ind,ng+1); /* +1 for possible rest group */
2157     for(i=0; (i<ng); i++)
2158     {
2159         /* Lookup the group name in the block structure */
2160         gid = search_string(ptrs[i],block->nr,gnames);
2161         if ((grptp != egrptpONE) || (i == 0))
2162         {
2163             grps->nm_ind[grps->nr++]=gid;
2164         }
2165         if (debug)
2166         {
2167             fprintf(debug,"Found gid %d for group %s\n",gid,ptrs[i]);
2168         }
2169
2170         /* Now go over the atoms in the group */
2171         for(j=block->index[gid]; (j<block->index[gid+1]); j++)
2172         {
2173
2174             aj=block->a[j];
2175
2176             /* Range checking */
2177             if ((aj < 0) || (aj >= natoms))
2178             {
2179                 gmx_fatal(FARGS,"Invalid atom number %d in indexfile",aj);
2180             }
2181             /* Lookup up the old group number */
2182             ognr = cbuf[aj];
2183             if (ognr != NOGID)
2184             {
2185                 gmx_fatal(FARGS,"Atom %d in multiple %s groups (%d and %d)",
2186                           aj+1,title,ognr+1,i+1);
2187             }
2188             else
2189             {
2190                 /* Store the group number in buffer */
2191                 if (grptp == egrptpONE)
2192                 {
2193                     cbuf[aj] = 0;
2194                 }
2195                 else
2196                 {
2197                     cbuf[aj] = i;
2198                 }
2199                 ntot++;
2200             }
2201         }
2202     }
2203
2204     /* Now check whether we have done all atoms */
2205     bRest = FALSE;
2206     if (ntot != natoms)
2207     {
2208         if (grptp == egrptpALL)
2209         {
2210             gmx_fatal(FARGS,"%d atoms are not part of any of the %s groups",
2211                       natoms-ntot,title);
2212         }
2213         else if (grptp == egrptpPART)
2214         {
2215             sprintf(warn_buf,"%d atoms are not part of any of the %s groups",
2216                     natoms-ntot,title);
2217             warning_note(wi,warn_buf);
2218         }
2219         /* Assign all atoms currently unassigned to a rest group */
2220         for(j=0; (j<natoms); j++)
2221         {
2222             if (cbuf[j] == NOGID)
2223             {
2224                 cbuf[j] = grps->nr;
2225                 bRest = TRUE;
2226             }
2227         }
2228         if (grptp != egrptpPART)
2229         {
2230             if (bVerbose)
2231             {
2232                 fprintf(stderr,
2233                         "Making dummy/rest group for %s containing %d elements\n",
2234                         title,natoms-ntot);
2235             }
2236             /* Add group name "rest" */
2237             grps->nm_ind[grps->nr] = restnm;
2238
2239             /* Assign the rest name to all atoms not currently assigned to a group */
2240             for(j=0; (j<natoms); j++)
2241             {
2242                 if (cbuf[j] == NOGID)
2243                 {
2244                     cbuf[j] = grps->nr;
2245                 }
2246             }
2247             grps->nr++;
2248         }
2249     }
2250
2251     if (grps->nr == 1 && (ntot == 0 || ntot == natoms))
2252     {
2253         /* All atoms are part of one (or no) group, no index required */
2254         groups->ngrpnr[gtype] = 0;
2255         groups->grpnr[gtype]  = NULL;
2256     }
2257     else
2258     {
2259         groups->ngrpnr[gtype] = natoms;
2260         snew(groups->grpnr[gtype],natoms);
2261         for(j=0; (j<natoms); j++)
2262         {
2263             groups->grpnr[gtype][j] = cbuf[j];
2264         }
2265     }
2266
2267     sfree(cbuf);
2268
2269     return (bRest && grptp == egrptpPART);
2270 }
2271
2272 static void calc_nrdf(gmx_mtop_t *mtop,t_inputrec *ir,char **gnames)
2273 {
2274   t_grpopts *opts;
2275   gmx_groups_t *groups;
2276   t_pull  *pull;
2277   int     natoms,ai,aj,i,j,d,g,imin,jmin,nc;
2278   t_iatom *ia;
2279   int     *nrdf2,*na_vcm,na_tot;
2280   double  *nrdf_tc,*nrdf_vcm,nrdf_uc,n_sub=0;
2281   gmx_mtop_atomloop_all_t aloop;
2282   t_atom  *atom;
2283   int     mb,mol,ftype,as;
2284   gmx_molblock_t *molb;
2285   gmx_moltype_t *molt;
2286
2287   /* Calculate nrdf.
2288    * First calc 3xnr-atoms for each group
2289    * then subtract half a degree of freedom for each constraint
2290    *
2291    * Only atoms and nuclei contribute to the degrees of freedom...
2292    */
2293
2294   opts = &ir->opts;
2295
2296   groups = &mtop->groups;
2297   natoms = mtop->natoms;
2298
2299   /* Allocate one more for a possible rest group */
2300   /* We need to sum degrees of freedom into doubles,
2301    * since floats give too low nrdf's above 3 million atoms.
2302    */
2303   snew(nrdf_tc,groups->grps[egcTC].nr+1);
2304   snew(nrdf_vcm,groups->grps[egcVCM].nr+1);
2305   snew(na_vcm,groups->grps[egcVCM].nr+1);
2306
2307   for(i=0; i<groups->grps[egcTC].nr; i++)
2308     nrdf_tc[i] = 0;
2309   for(i=0; i<groups->grps[egcVCM].nr+1; i++)
2310     nrdf_vcm[i] = 0;
2311
2312   snew(nrdf2,natoms);
2313   aloop = gmx_mtop_atomloop_all_init(mtop);
2314   while (gmx_mtop_atomloop_all_next(aloop,&i,&atom)) {
2315     nrdf2[i] = 0;
2316     if (atom->ptype == eptAtom || atom->ptype == eptNucleus) {
2317       g = ggrpnr(groups,egcFREEZE,i);
2318       /* Double count nrdf for particle i */
2319       for(d=0; d<DIM; d++) {
2320         if (opts->nFreeze[g][d] == 0) {
2321           nrdf2[i] += 2;
2322         }
2323       }
2324       nrdf_tc [ggrpnr(groups,egcTC ,i)] += 0.5*nrdf2[i];
2325       nrdf_vcm[ggrpnr(groups,egcVCM,i)] += 0.5*nrdf2[i];
2326     }
2327   }
2328
2329   as = 0;
2330   for(mb=0; mb<mtop->nmolblock; mb++) {
2331     molb = &mtop->molblock[mb];
2332     molt = &mtop->moltype[molb->type];
2333     atom = molt->atoms.atom;
2334     for(mol=0; mol<molb->nmol; mol++) {
2335       for (ftype=F_CONSTR; ftype<=F_CONSTRNC; ftype++) {
2336         ia = molt->ilist[ftype].iatoms;
2337         for(i=0; i<molt->ilist[ftype].nr; ) {
2338           /* Subtract degrees of freedom for the constraints,
2339            * if the particles still have degrees of freedom left.
2340            * If one of the particles is a vsite or a shell, then all
2341            * constraint motion will go there, but since they do not
2342            * contribute to the constraints the degrees of freedom do not
2343            * change.
2344            */
2345           ai = as + ia[1];
2346           aj = as + ia[2];
2347           if (((atom[ia[1]].ptype == eptNucleus) ||
2348                (atom[ia[1]].ptype == eptAtom)) &&
2349               ((atom[ia[2]].ptype == eptNucleus) ||
2350                (atom[ia[2]].ptype == eptAtom))) {
2351             if (nrdf2[ai] > 0)
2352               jmin = 1;
2353             else
2354               jmin = 2;
2355             if (nrdf2[aj] > 0)
2356               imin = 1;
2357             else
2358               imin = 2;
2359             imin = min(imin,nrdf2[ai]);
2360             jmin = min(jmin,nrdf2[aj]);
2361             nrdf2[ai] -= imin;
2362             nrdf2[aj] -= jmin;
2363             nrdf_tc [ggrpnr(groups,egcTC ,ai)] -= 0.5*imin;
2364             nrdf_tc [ggrpnr(groups,egcTC ,aj)] -= 0.5*jmin;
2365             nrdf_vcm[ggrpnr(groups,egcVCM,ai)] -= 0.5*imin;
2366             nrdf_vcm[ggrpnr(groups,egcVCM,aj)] -= 0.5*jmin;
2367           }
2368           ia += interaction_function[ftype].nratoms+1;
2369           i  += interaction_function[ftype].nratoms+1;
2370         }
2371       }
2372       ia = molt->ilist[F_SETTLE].iatoms;
2373       for(i=0; i<molt->ilist[F_SETTLE].nr; ) {
2374         /* Subtract 1 dof from every atom in the SETTLE */
2375         for(j=0; j<3; j++) {
2376       ai = as + ia[1+j];
2377           imin = min(2,nrdf2[ai]);
2378           nrdf2[ai] -= imin;
2379           nrdf_tc [ggrpnr(groups,egcTC ,ai)] -= 0.5*imin;
2380           nrdf_vcm[ggrpnr(groups,egcVCM,ai)] -= 0.5*imin;
2381         }
2382         ia += 4;
2383         i  += 4;
2384       }
2385       as += molt->atoms.nr;
2386     }
2387   }
2388
2389   if (ir->ePull == epullCONSTRAINT) {
2390     /* Correct nrdf for the COM constraints.
2391      * We correct using the TC and VCM group of the first atom
2392      * in the reference and pull group. If atoms in one pull group
2393      * belong to different TC or VCM groups it is anyhow difficult
2394      * to determine the optimal nrdf assignment.
2395      */
2396     pull = ir->pull;
2397     if (pull->eGeom == epullgPOS) {
2398       nc = 0;
2399       for(i=0; i<DIM; i++) {
2400         if (pull->dim[i])
2401           nc++;
2402       }
2403     } else {
2404       nc = 1;
2405     }
2406     for(i=0; i<pull->ngrp; i++) {
2407       imin = 2*nc;
2408       if (pull->grp[0].nat > 0) {
2409         /* Subtract 1/2 dof from the reference group */
2410         ai = pull->grp[0].ind[0];
2411         if (nrdf_tc[ggrpnr(groups,egcTC,ai)] > 1) {
2412           nrdf_tc [ggrpnr(groups,egcTC ,ai)] -= 0.5;
2413           nrdf_vcm[ggrpnr(groups,egcVCM,ai)] -= 0.5;
2414           imin--;
2415         }
2416       }
2417       /* Subtract 1/2 dof from the pulled group */
2418       ai = pull->grp[1+i].ind[0];
2419       nrdf_tc [ggrpnr(groups,egcTC ,ai)] -= 0.5*imin;
2420       nrdf_vcm[ggrpnr(groups,egcVCM,ai)] -= 0.5*imin;
2421       if (nrdf_tc[ggrpnr(groups,egcTC,ai)] < 0)
2422         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)]]);
2423     }
2424   }
2425
2426   if (ir->nstcomm != 0) {
2427     /* Subtract 3 from the number of degrees of freedom in each vcm group
2428      * when com translation is removed and 6 when rotation is removed
2429      * as well.
2430      */
2431     switch (ir->comm_mode) {
2432     case ecmLINEAR:
2433       n_sub = ndof_com(ir);
2434       break;
2435     case ecmANGULAR:
2436       n_sub = 6;
2437       break;
2438     default:
2439       n_sub = 0;
2440       gmx_incons("Checking comm_mode");
2441     }
2442
2443     for(i=0; i<groups->grps[egcTC].nr; i++) {
2444       /* Count the number of atoms of TC group i for every VCM group */
2445       for(j=0; j<groups->grps[egcVCM].nr+1; j++)
2446         na_vcm[j] = 0;
2447       na_tot = 0;
2448       for(ai=0; ai<natoms; ai++)
2449         if (ggrpnr(groups,egcTC,ai) == i) {
2450           na_vcm[ggrpnr(groups,egcVCM,ai)]++;
2451           na_tot++;
2452         }
2453       /* Correct for VCM removal according to the fraction of each VCM
2454        * group present in this TC group.
2455        */
2456       nrdf_uc = nrdf_tc[i];
2457       if (debug) {
2458         fprintf(debug,"T-group[%d] nrdf_uc = %g, n_sub = %g\n",
2459                 i,nrdf_uc,n_sub);
2460       }
2461       nrdf_tc[i] = 0;
2462       for(j=0; j<groups->grps[egcVCM].nr+1; j++) {
2463         if (nrdf_vcm[j] > n_sub) {
2464           nrdf_tc[i] += nrdf_uc*((double)na_vcm[j]/(double)na_tot)*
2465             (nrdf_vcm[j] - n_sub)/nrdf_vcm[j];
2466         }
2467         if (debug) {
2468           fprintf(debug,"  nrdf_vcm[%d] = %g, nrdf = %g\n",
2469                   j,nrdf_vcm[j],nrdf_tc[i]);
2470         }
2471       }
2472     }
2473   }
2474   for(i=0; (i<groups->grps[egcTC].nr); i++) {
2475     opts->nrdf[i] = nrdf_tc[i];
2476     if (opts->nrdf[i] < 0)
2477       opts->nrdf[i] = 0;
2478     fprintf(stderr,
2479             "Number of degrees of freedom in T-Coupling group %s is %.2f\n",
2480             gnames[groups->grps[egcTC].nm_ind[i]],opts->nrdf[i]);
2481   }
2482
2483   sfree(nrdf2);
2484   sfree(nrdf_tc);
2485   sfree(nrdf_vcm);
2486   sfree(na_vcm);
2487 }
2488
2489 static void decode_cos(char *s,t_cosines *cosine,gmx_bool bTime)
2490 {
2491   char   *t;
2492   char   format[STRLEN],f1[STRLEN];
2493   double a,phi;
2494   int    i;
2495
2496   t=strdup(s);
2497   trim(t);
2498
2499   cosine->n=0;
2500   cosine->a=NULL;
2501   cosine->phi=NULL;
2502   if (strlen(t)) {
2503     sscanf(t,"%d",&(cosine->n));
2504     if (cosine->n <= 0) {
2505       cosine->n=0;
2506     } else {
2507       snew(cosine->a,cosine->n);
2508       snew(cosine->phi,cosine->n);
2509
2510       sprintf(format,"%%*d");
2511       for(i=0; (i<cosine->n); i++) {
2512         strcpy(f1,format);
2513         strcat(f1,"%lf%lf");
2514         if (sscanf(t,f1,&a,&phi) < 2)
2515           gmx_fatal(FARGS,"Invalid input for electric field shift: '%s'",t);
2516         cosine->a[i]=a;
2517         cosine->phi[i]=phi;
2518         strcat(format,"%*lf%*lf");
2519       }
2520     }
2521   }
2522   sfree(t);
2523 }
2524
2525 static gmx_bool do_egp_flag(t_inputrec *ir,gmx_groups_t *groups,
2526                         const char *option,const char *val,int flag)
2527 {
2528   /* The maximum number of energy group pairs would be MAXPTR*(MAXPTR+1)/2.
2529    * But since this is much larger than STRLEN, such a line can not be parsed.
2530    * The real maximum is the number of names that fit in a string: STRLEN/2.
2531    */
2532 #define EGP_MAX (STRLEN/2)
2533   int  nelem,i,j,k,nr;
2534   char *names[EGP_MAX];
2535   char ***gnames;
2536   gmx_bool bSet;
2537
2538   gnames = groups->grpname;
2539
2540   nelem = str_nelem(val,EGP_MAX,names);
2541   if (nelem % 2 != 0)
2542     gmx_fatal(FARGS,"The number of groups for %s is odd",option);
2543   nr = groups->grps[egcENER].nr;
2544   bSet = FALSE;
2545   for(i=0; i<nelem/2; i++) {
2546     j = 0;
2547     while ((j < nr) &&
2548            gmx_strcasecmp(names[2*i],*(gnames[groups->grps[egcENER].nm_ind[j]])))
2549       j++;
2550     if (j == nr)
2551       gmx_fatal(FARGS,"%s in %s is not an energy group\n",
2552                   names[2*i],option);
2553     k = 0;
2554     while ((k < nr) &&
2555            gmx_strcasecmp(names[2*i+1],*(gnames[groups->grps[egcENER].nm_ind[k]])))
2556       k++;
2557     if (k==nr)
2558       gmx_fatal(FARGS,"%s in %s is not an energy group\n",
2559               names[2*i+1],option);
2560     if ((j < nr) && (k < nr)) {
2561       ir->opts.egp_flags[nr*j+k] |= flag;
2562       ir->opts.egp_flags[nr*k+j] |= flag;
2563       bSet = TRUE;
2564     }
2565   }
2566
2567   return bSet;
2568 }
2569
2570 void do_index(const char* mdparin, const char *ndx,
2571               gmx_mtop_t *mtop,
2572               gmx_bool bVerbose,
2573               t_inputrec *ir,rvec *v,
2574               warninp_t wi)
2575 {
2576   t_blocka *grps;
2577   gmx_groups_t *groups;
2578   int     natoms;
2579   t_symtab *symtab;
2580   t_atoms atoms_all;
2581   char    warnbuf[STRLEN],**gnames;
2582   int     nr,ntcg,ntau_t,nref_t,nacc,nofg,nSA,nSA_points,nSA_time,nSA_temp;
2583   real    tau_min;
2584   int     nstcmin;
2585   int     nacg,nfreeze,nfrdim,nenergy,nvcm,nuser;
2586   char    *ptr1[MAXPTR],*ptr2[MAXPTR],*ptr3[MAXPTR];
2587   int     i,j,k,restnm;
2588   real    SAtime;
2589   gmx_bool    bExcl,bTable,bSetTCpar,bAnneal,bRest;
2590   int     nQMmethod,nQMbasis,nQMcharge,nQMmult,nbSH,nCASorb,nCASelec,
2591     nSAon,nSAoff,nSAsteps,nQMg,nbOPT,nbTS;
2592   char    warn_buf[STRLEN];
2593
2594   if (bVerbose)
2595     fprintf(stderr,"processing index file...\n");
2596   debug_gmx();
2597   if (ndx == NULL) {
2598     snew(grps,1);
2599     snew(grps->index,1);
2600     snew(gnames,1);
2601     atoms_all = gmx_mtop_global_atoms(mtop);
2602     analyse(&atoms_all,grps,&gnames,FALSE,TRUE);
2603     free_t_atoms(&atoms_all,FALSE);
2604   } else {
2605     grps = init_index(ndx,&gnames);
2606   }
2607
2608   groups = &mtop->groups;
2609   natoms = mtop->natoms;
2610   symtab = &mtop->symtab;
2611
2612   snew(groups->grpname,grps->nr+1);
2613
2614   for(i=0; (i<grps->nr); i++) {
2615     groups->grpname[i] = put_symtab(symtab,gnames[i]);
2616   }
2617   groups->grpname[i] = put_symtab(symtab,"rest");
2618   restnm=i;
2619   srenew(gnames,grps->nr+1);
2620   gnames[restnm] = *(groups->grpname[i]);
2621   groups->ngrpname = grps->nr+1;
2622
2623   set_warning_line(wi,mdparin,-1);
2624
2625   ntau_t = str_nelem(tau_t,MAXPTR,ptr1);
2626   nref_t = str_nelem(ref_t,MAXPTR,ptr2);
2627   ntcg   = str_nelem(tcgrps,MAXPTR,ptr3);
2628   if ((ntau_t != ntcg) || (nref_t != ntcg)) {
2629     gmx_fatal(FARGS,"Invalid T coupling input: %d groups, %d ref-t values and "
2630                 "%d tau-t values",ntcg,nref_t,ntau_t);
2631   }
2632
2633   bSetTCpar = (ir->etc || EI_SD(ir->eI) || ir->eI==eiBD || EI_TPI(ir->eI));
2634   do_numbering(natoms,groups,ntcg,ptr3,grps,gnames,egcTC,
2635                restnm,bSetTCpar ? egrptpALL : egrptpALL_GENREST,bVerbose,wi);
2636   nr = groups->grps[egcTC].nr;
2637   ir->opts.ngtc = nr;
2638   snew(ir->opts.nrdf,nr);
2639   snew(ir->opts.tau_t,nr);
2640   snew(ir->opts.ref_t,nr);
2641   if (ir->eI==eiBD && ir->bd_fric==0) {
2642     fprintf(stderr,"bd-fric=0, so tau-t will be used as the inverse friction constant(s)\n");
2643   }
2644
2645   if (bSetTCpar)
2646   {
2647       if (nr != nref_t)
2648       {
2649           gmx_fatal(FARGS,"Not enough ref-t and tau-t values!");
2650       }
2651
2652       tau_min = 1e20;
2653       for(i=0; (i<nr); i++)
2654       {
2655           ir->opts.tau_t[i] = strtod(ptr1[i],NULL);
2656           if ((ir->eI == eiBD || ir->eI == eiSD2) && ir->opts.tau_t[i] <= 0)
2657           {
2658               sprintf(warn_buf,"With integrator %s tau-t should be larger than 0",ei_names[ir->eI]);
2659               warning_error(wi,warn_buf);
2660           }
2661
2662           if (ir->etc != etcVRESCALE && ir->opts.tau_t[i] == 0)
2663           {
2664               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.");
2665           }
2666
2667           if (ir->opts.tau_t[i] >= 0)
2668           {
2669               tau_min = min(tau_min,ir->opts.tau_t[i]);
2670           }
2671       }
2672       if (ir->etc != etcNO && ir->nsttcouple == -1)
2673       {
2674             ir->nsttcouple = ir_optimal_nsttcouple(ir);
2675       }
2676
2677       if (EI_VV(ir->eI))
2678       {
2679           if ((ir->etc==etcNOSEHOOVER) && (ir->epc==epcBERENDSEN)) {
2680               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");
2681           }
2682           if ((ir->epc==epcMTTK) && (ir->etc>etcNO))
2683           {
2684               int mincouple;
2685               mincouple = ir->nsttcouple;
2686               if (ir->nstpcouple < mincouple)
2687               {
2688                   mincouple = ir->nstpcouple;
2689               }
2690               ir->nstpcouple = mincouple;
2691               ir->nsttcouple = mincouple;
2692               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);
2693               warning_note(wi,warn_buf);
2694           }
2695       }
2696       /* velocity verlet with averaged kinetic energy KE = 0.5*(v(t+1/2) - v(t-1/2)) is implemented
2697          primarily for testing purposes, and does not work with temperature coupling other than 1 */
2698
2699       if (ETC_ANDERSEN(ir->etc)) {
2700           if (ir->nsttcouple != 1) {
2701               ir->nsttcouple = 1;
2702               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");
2703               warning_note(wi,warn_buf);
2704           }
2705       }
2706       nstcmin = tcouple_min_integration_steps(ir->etc);
2707       if (nstcmin > 1)
2708       {
2709           if (tau_min/(ir->delta_t*ir->nsttcouple) < nstcmin)
2710           {
2711               sprintf(warn_buf,"For proper integration of the %s thermostat, tau-t (%g) should be at least %d times larger than nsttcouple*dt (%g)",
2712                       ETCOUPLTYPE(ir->etc),
2713                       tau_min,nstcmin,
2714                       ir->nsttcouple*ir->delta_t);
2715               warning(wi,warn_buf);
2716           }
2717       }
2718       for(i=0; (i<nr); i++)
2719       {
2720           ir->opts.ref_t[i] = strtod(ptr2[i],NULL);
2721           if (ir->opts.ref_t[i] < 0)
2722           {
2723               gmx_fatal(FARGS,"ref-t for group %d negative",i);
2724           }
2725       }
2726       /* set the lambda mc temperature to the md integrator temperature (which should be defined
2727          if we are in this conditional) if mc_temp is negative */
2728       if (ir->expandedvals->mc_temp < 0)
2729       {
2730           ir->expandedvals->mc_temp = ir->opts.ref_t[0];  /*for now, set to the first reft */
2731       }
2732   }
2733
2734   /* Simulated annealing for each group. There are nr groups */
2735   nSA = str_nelem(anneal,MAXPTR,ptr1);
2736   if (nSA == 1 && (ptr1[0][0]=='n' || ptr1[0][0]=='N'))
2737      nSA = 0;
2738   if(nSA>0 && nSA != nr)
2739     gmx_fatal(FARGS,"Not enough annealing values: %d (for %d groups)\n",nSA,nr);
2740   else {
2741     snew(ir->opts.annealing,nr);
2742     snew(ir->opts.anneal_npoints,nr);
2743     snew(ir->opts.anneal_time,nr);
2744     snew(ir->opts.anneal_temp,nr);
2745     for(i=0;i<nr;i++) {
2746       ir->opts.annealing[i]=eannNO;
2747       ir->opts.anneal_npoints[i]=0;
2748       ir->opts.anneal_time[i]=NULL;
2749       ir->opts.anneal_temp[i]=NULL;
2750     }
2751     if (nSA > 0) {
2752       bAnneal=FALSE;
2753       for(i=0;i<nr;i++) {
2754         if(ptr1[i][0]=='n' || ptr1[i][0]=='N') {
2755           ir->opts.annealing[i]=eannNO;
2756         } else if(ptr1[i][0]=='s'|| ptr1[i][0]=='S') {
2757           ir->opts.annealing[i]=eannSINGLE;
2758           bAnneal=TRUE;
2759         } else if(ptr1[i][0]=='p'|| ptr1[i][0]=='P') {
2760           ir->opts.annealing[i]=eannPERIODIC;
2761           bAnneal=TRUE;
2762         }
2763       }
2764       if(bAnneal) {
2765         /* Read the other fields too */
2766         nSA_points = str_nelem(anneal_npoints,MAXPTR,ptr1);
2767         if(nSA_points!=nSA)
2768           gmx_fatal(FARGS,"Found %d annealing-npoints values for %d groups\n",nSA_points,nSA);
2769         for(k=0,i=0;i<nr;i++) {
2770           ir->opts.anneal_npoints[i]=strtol(ptr1[i],NULL,10);
2771           if(ir->opts.anneal_npoints[i]==1)
2772             gmx_fatal(FARGS,"Please specify at least a start and an end point for annealing\n");
2773           snew(ir->opts.anneal_time[i],ir->opts.anneal_npoints[i]);
2774           snew(ir->opts.anneal_temp[i],ir->opts.anneal_npoints[i]);
2775           k += ir->opts.anneal_npoints[i];
2776         }
2777
2778         nSA_time = str_nelem(anneal_time,MAXPTR,ptr1);
2779         if(nSA_time!=k)
2780           gmx_fatal(FARGS,"Found %d annealing-time values, wanter %d\n",nSA_time,k);
2781         nSA_temp = str_nelem(anneal_temp,MAXPTR,ptr2);
2782         if(nSA_temp!=k)
2783           gmx_fatal(FARGS,"Found %d annealing-temp values, wanted %d\n",nSA_temp,k);
2784
2785         for(i=0,k=0;i<nr;i++) {
2786
2787           for(j=0;j<ir->opts.anneal_npoints[i];j++) {
2788             ir->opts.anneal_time[i][j]=strtod(ptr1[k],NULL);
2789             ir->opts.anneal_temp[i][j]=strtod(ptr2[k],NULL);
2790             if(j==0) {
2791               if(ir->opts.anneal_time[i][0] > (ir->init_t+GMX_REAL_EPS))
2792                 gmx_fatal(FARGS,"First time point for annealing > init_t.\n");
2793             } else {
2794               /* j>0 */
2795               if(ir->opts.anneal_time[i][j]<ir->opts.anneal_time[i][j-1])
2796                 gmx_fatal(FARGS,"Annealing timepoints out of order: t=%f comes after t=%f\n",
2797                             ir->opts.anneal_time[i][j],ir->opts.anneal_time[i][j-1]);
2798             }
2799             if(ir->opts.anneal_temp[i][j]<0)
2800               gmx_fatal(FARGS,"Found negative temperature in annealing: %f\n",ir->opts.anneal_temp[i][j]);
2801             k++;
2802           }
2803         }
2804         /* Print out some summary information, to make sure we got it right */
2805         for(i=0,k=0;i<nr;i++) {
2806           if(ir->opts.annealing[i]!=eannNO) {
2807             j = groups->grps[egcTC].nm_ind[i];
2808             fprintf(stderr,"Simulated annealing for group %s: %s, %d timepoints\n",
2809                     *(groups->grpname[j]),eann_names[ir->opts.annealing[i]],
2810                     ir->opts.anneal_npoints[i]);
2811             fprintf(stderr,"Time (ps)   Temperature (K)\n");
2812             /* All terms except the last one */
2813             for(j=0;j<(ir->opts.anneal_npoints[i]-1);j++)
2814                 fprintf(stderr,"%9.1f      %5.1f\n",ir->opts.anneal_time[i][j],ir->opts.anneal_temp[i][j]);
2815
2816             /* Finally the last one */
2817             j = ir->opts.anneal_npoints[i]-1;
2818             if(ir->opts.annealing[i]==eannSINGLE)
2819               fprintf(stderr,"%9.1f-     %5.1f\n",ir->opts.anneal_time[i][j],ir->opts.anneal_temp[i][j]);
2820             else {
2821               fprintf(stderr,"%9.1f      %5.1f\n",ir->opts.anneal_time[i][j],ir->opts.anneal_temp[i][j]);
2822               if(fabs(ir->opts.anneal_temp[i][j]-ir->opts.anneal_temp[i][0])>GMX_REAL_EPS)
2823                 warning_note(wi,"There is a temperature jump when your annealing loops back.\n");
2824             }
2825           }
2826         }
2827       }
2828     }
2829   }
2830
2831   if (ir->ePull != epullNO) {
2832     make_pull_groups(ir->pull,pull_grp,grps,gnames);
2833   }
2834
2835   if (ir->bRot) {
2836     make_rotation_groups(ir->rot,rot_grp,grps,gnames);
2837   }
2838
2839   nacc = str_nelem(acc,MAXPTR,ptr1);
2840   nacg = str_nelem(accgrps,MAXPTR,ptr2);
2841   if (nacg*DIM != nacc)
2842     gmx_fatal(FARGS,"Invalid Acceleration input: %d groups and %d acc. values",
2843                 nacg,nacc);
2844   do_numbering(natoms,groups,nacg,ptr2,grps,gnames,egcACC,
2845                restnm,egrptpALL_GENREST,bVerbose,wi);
2846   nr = groups->grps[egcACC].nr;
2847   snew(ir->opts.acc,nr);
2848   ir->opts.ngacc=nr;
2849
2850   for(i=k=0; (i<nacg); i++)
2851     for(j=0; (j<DIM); j++,k++)
2852       ir->opts.acc[i][j]=strtod(ptr1[k],NULL);
2853   for( ;(i<nr); i++)
2854     for(j=0; (j<DIM); j++)
2855       ir->opts.acc[i][j]=0;
2856
2857   nfrdim  = str_nelem(frdim,MAXPTR,ptr1);
2858   nfreeze = str_nelem(freeze,MAXPTR,ptr2);
2859   if (nfrdim != DIM*nfreeze)
2860     gmx_fatal(FARGS,"Invalid Freezing input: %d groups and %d freeze values",
2861                 nfreeze,nfrdim);
2862   do_numbering(natoms,groups,nfreeze,ptr2,grps,gnames,egcFREEZE,
2863                restnm,egrptpALL_GENREST,bVerbose,wi);
2864   nr = groups->grps[egcFREEZE].nr;
2865   ir->opts.ngfrz=nr;
2866   snew(ir->opts.nFreeze,nr);
2867   for(i=k=0; (i<nfreeze); i++)
2868     for(j=0; (j<DIM); j++,k++) {
2869       ir->opts.nFreeze[i][j]=(gmx_strncasecmp(ptr1[k],"Y",1)==0);
2870       if (!ir->opts.nFreeze[i][j]) {
2871         if (gmx_strncasecmp(ptr1[k],"N",1) != 0) {
2872           sprintf(warnbuf,"Please use Y(ES) or N(O) for freezedim only "
2873                   "(not %s)", ptr1[k]);
2874           warning(wi,warn_buf);
2875         }
2876       }
2877     }
2878   for( ; (i<nr); i++)
2879     for(j=0; (j<DIM); j++)
2880       ir->opts.nFreeze[i][j]=0;
2881
2882   nenergy=str_nelem(energy,MAXPTR,ptr1);
2883   do_numbering(natoms,groups,nenergy,ptr1,grps,gnames,egcENER,
2884                restnm,egrptpALL_GENREST,bVerbose,wi);
2885   add_wall_energrps(groups,ir->nwall,symtab);
2886   ir->opts.ngener = groups->grps[egcENER].nr;
2887   nvcm=str_nelem(vcm,MAXPTR,ptr1);
2888   bRest =
2889     do_numbering(natoms,groups,nvcm,ptr1,grps,gnames,egcVCM,
2890                  restnm,nvcm==0 ? egrptpALL_GENREST : egrptpPART,bVerbose,wi);
2891   if (bRest) {
2892     warning(wi,"Some atoms are not part of any center of mass motion removal group.\n"
2893             "This may lead to artifacts.\n"
2894             "In most cases one should use one group for the whole system.");
2895   }
2896
2897   /* Now we have filled the freeze struct, so we can calculate NRDF */
2898   calc_nrdf(mtop,ir,gnames);
2899
2900   if (v && NULL) {
2901     real fac,ntot=0;
2902
2903     /* Must check per group! */
2904     for(i=0; (i<ir->opts.ngtc); i++)
2905       ntot += ir->opts.nrdf[i];
2906     if (ntot != (DIM*natoms)) {
2907       fac = sqrt(ntot/(DIM*natoms));
2908       if (bVerbose)
2909         fprintf(stderr,"Scaling velocities by a factor of %.3f to account for constraints\n"
2910                 "and removal of center of mass motion\n",fac);
2911       for(i=0; (i<natoms); i++)
2912         svmul(fac,v[i],v[i]);
2913     }
2914   }
2915
2916   nuser=str_nelem(user1,MAXPTR,ptr1);
2917   do_numbering(natoms,groups,nuser,ptr1,grps,gnames,egcUser1,
2918                restnm,egrptpALL_GENREST,bVerbose,wi);
2919   nuser=str_nelem(user2,MAXPTR,ptr1);
2920   do_numbering(natoms,groups,nuser,ptr1,grps,gnames,egcUser2,
2921                restnm,egrptpALL_GENREST,bVerbose,wi);
2922   nuser=str_nelem(xtc_grps,MAXPTR,ptr1);
2923   do_numbering(natoms,groups,nuser,ptr1,grps,gnames,egcXTC,
2924                restnm,egrptpONE,bVerbose,wi);
2925   nofg = str_nelem(orirefitgrp,MAXPTR,ptr1);
2926   do_numbering(natoms,groups,nofg,ptr1,grps,gnames,egcORFIT,
2927                restnm,egrptpALL_GENREST,bVerbose,wi);
2928
2929   /* QMMM input processing */
2930   nQMg          = str_nelem(QMMM,MAXPTR,ptr1);
2931   nQMmethod     = str_nelem(QMmethod,MAXPTR,ptr2);
2932   nQMbasis      = str_nelem(QMbasis,MAXPTR,ptr3);
2933   if((nQMmethod != nQMg)||(nQMbasis != nQMg)){
2934     gmx_fatal(FARGS,"Invalid QMMM input: %d groups %d basissets"
2935               " and %d methods\n",nQMg,nQMbasis,nQMmethod);
2936   }
2937   /* group rest, if any, is always MM! */
2938   do_numbering(natoms,groups,nQMg,ptr1,grps,gnames,egcQMMM,
2939                restnm,egrptpALL_GENREST,bVerbose,wi);
2940   nr = nQMg; /*atoms->grps[egcQMMM].nr;*/
2941   ir->opts.ngQM = nQMg;
2942   snew(ir->opts.QMmethod,nr);
2943   snew(ir->opts.QMbasis,nr);
2944   for(i=0;i<nr;i++){
2945     /* input consists of strings: RHF CASSCF PM3 .. These need to be
2946      * converted to the corresponding enum in names.c
2947      */
2948     ir->opts.QMmethod[i] = search_QMstring(ptr2[i],eQMmethodNR,
2949                                            eQMmethod_names);
2950     ir->opts.QMbasis[i]  = search_QMstring(ptr3[i],eQMbasisNR,
2951                                            eQMbasis_names);
2952
2953   }
2954   nQMmult   = str_nelem(QMmult,MAXPTR,ptr1);
2955   nQMcharge = str_nelem(QMcharge,MAXPTR,ptr2);
2956   nbSH      = str_nelem(bSH,MAXPTR,ptr3);
2957   snew(ir->opts.QMmult,nr);
2958   snew(ir->opts.QMcharge,nr);
2959   snew(ir->opts.bSH,nr);
2960
2961   for(i=0;i<nr;i++){
2962     ir->opts.QMmult[i]   = strtol(ptr1[i],NULL,10);
2963     ir->opts.QMcharge[i] = strtol(ptr2[i],NULL,10);
2964     ir->opts.bSH[i]      = (gmx_strncasecmp(ptr3[i],"Y",1)==0);
2965   }
2966
2967   nCASelec  = str_nelem(CASelectrons,MAXPTR,ptr1);
2968   nCASorb   = str_nelem(CASorbitals,MAXPTR,ptr2);
2969   snew(ir->opts.CASelectrons,nr);
2970   snew(ir->opts.CASorbitals,nr);
2971   for(i=0;i<nr;i++){
2972     ir->opts.CASelectrons[i]= strtol(ptr1[i],NULL,10);
2973     ir->opts.CASorbitals[i] = strtol(ptr2[i],NULL,10);
2974   }
2975   /* special optimization options */
2976
2977   nbOPT = str_nelem(bOPT,MAXPTR,ptr1);
2978   nbTS = str_nelem(bTS,MAXPTR,ptr2);
2979   snew(ir->opts.bOPT,nr);
2980   snew(ir->opts.bTS,nr);
2981   for(i=0;i<nr;i++){
2982     ir->opts.bOPT[i] = (gmx_strncasecmp(ptr1[i],"Y",1)==0);
2983     ir->opts.bTS[i]  = (gmx_strncasecmp(ptr2[i],"Y",1)==0);
2984   }
2985   nSAon     = str_nelem(SAon,MAXPTR,ptr1);
2986   nSAoff    = str_nelem(SAoff,MAXPTR,ptr2);
2987   nSAsteps  = str_nelem(SAsteps,MAXPTR,ptr3);
2988   snew(ir->opts.SAon,nr);
2989   snew(ir->opts.SAoff,nr);
2990   snew(ir->opts.SAsteps,nr);
2991
2992   for(i=0;i<nr;i++){
2993     ir->opts.SAon[i]    = strtod(ptr1[i],NULL);
2994     ir->opts.SAoff[i]   = strtod(ptr2[i],NULL);
2995     ir->opts.SAsteps[i] = strtol(ptr3[i],NULL,10);
2996   }
2997   /* end of QMMM input */
2998
2999   if (bVerbose)
3000     for(i=0; (i<egcNR); i++) {
3001       fprintf(stderr,"%-16s has %d element(s):",gtypes[i],groups->grps[i].nr);
3002       for(j=0; (j<groups->grps[i].nr); j++)
3003         fprintf(stderr," %s",*(groups->grpname[groups->grps[i].nm_ind[j]]));
3004       fprintf(stderr,"\n");
3005     }
3006
3007   nr = groups->grps[egcENER].nr;
3008   snew(ir->opts.egp_flags,nr*nr);
3009
3010   bExcl = do_egp_flag(ir,groups,"energygrp-excl",egpexcl,EGP_EXCL);
3011     if (bExcl && ir->cutoff_scheme == ecutsVERLET)
3012     {
3013         warning_error(wi,"Energy group exclusions are not (yet) implemented for the Verlet scheme");
3014     }
3015   if (bExcl && EEL_FULL(ir->coulombtype))
3016     warning(wi,"Can not exclude the lattice Coulomb energy between energy groups");
3017
3018   bTable = do_egp_flag(ir,groups,"energygrp-table",egptable,EGP_TABLE);
3019   if (bTable && !(ir->vdwtype == evdwUSER) &&
3020       !(ir->coulombtype == eelUSER) && !(ir->coulombtype == eelPMEUSER) &&
3021       !(ir->coulombtype == eelPMEUSERSWITCH))
3022     gmx_fatal(FARGS,"Can only have energy group pair tables in combination with user tables for VdW and/or Coulomb");
3023
3024   decode_cos(efield_x,&(ir->ex[XX]),FALSE);
3025   decode_cos(efield_xt,&(ir->et[XX]),TRUE);
3026   decode_cos(efield_y,&(ir->ex[YY]),FALSE);
3027   decode_cos(efield_yt,&(ir->et[YY]),TRUE);
3028   decode_cos(efield_z,&(ir->ex[ZZ]),FALSE);
3029   decode_cos(efield_zt,&(ir->et[ZZ]),TRUE);
3030
3031   if (ir->bAdress)
3032     do_adress_index(ir->adress,groups,gnames,&(ir->opts),wi);
3033
3034   for(i=0; (i<grps->nr); i++)
3035     sfree(gnames[i]);
3036   sfree(gnames);
3037   done_blocka(grps);
3038   sfree(grps);
3039
3040 }
3041
3042
3043
3044 static void check_disre(gmx_mtop_t *mtop)
3045 {
3046   gmx_ffparams_t *ffparams;
3047   t_functype *functype;
3048   t_iparams  *ip;
3049   int i,ndouble,ftype;
3050   int label,old_label;
3051
3052   if (gmx_mtop_ftype_count(mtop,F_DISRES) > 0) {
3053     ffparams  = &mtop->ffparams;
3054     functype  = ffparams->functype;
3055     ip        = ffparams->iparams;
3056     ndouble   = 0;
3057     old_label = -1;
3058     for(i=0; i<ffparams->ntypes; i++) {
3059       ftype = functype[i];
3060       if (ftype == F_DISRES) {
3061         label = ip[i].disres.label;
3062         if (label == old_label) {
3063           fprintf(stderr,"Distance restraint index %d occurs twice\n",label);
3064           ndouble++;
3065         }
3066         old_label = label;
3067       }
3068     }
3069     if (ndouble>0)
3070       gmx_fatal(FARGS,"Found %d double distance restraint indices,\n"
3071                 "probably the parameters for multiple pairs in one restraint "
3072                 "are not identical\n",ndouble);
3073   }
3074 }
3075
3076 static gmx_bool absolute_reference(t_inputrec *ir,gmx_mtop_t *sys,
3077                                    gmx_bool posres_only,
3078                                    ivec AbsRef)
3079 {
3080     int d,g,i;
3081     gmx_mtop_ilistloop_t iloop;
3082     t_ilist *ilist;
3083     int nmol;
3084     t_iparams *pr;
3085
3086     clear_ivec(AbsRef);
3087
3088     if (!posres_only)
3089     {
3090         /* Check the COM */
3091         for(d=0; d<DIM; d++)
3092         {
3093             AbsRef[d] = (d < ndof_com(ir) ? 0 : 1);
3094         }
3095         /* Check for freeze groups */
3096         for(g=0; g<ir->opts.ngfrz; g++)
3097         {
3098             for(d=0; d<DIM; d++)
3099             {
3100                 if (ir->opts.nFreeze[g][d] != 0)
3101                 {
3102                     AbsRef[d] = 1;
3103                 }
3104             }
3105         }
3106     }
3107
3108     /* Check for position restraints */
3109     iloop = gmx_mtop_ilistloop_init(sys);
3110     while (gmx_mtop_ilistloop_next(iloop,&ilist,&nmol))
3111     {
3112         if (nmol > 0 &&
3113             (AbsRef[XX] == 0 || AbsRef[YY] == 0 || AbsRef[ZZ] == 0))
3114         {
3115             for(i=0; i<ilist[F_POSRES].nr; i+=2)
3116             {
3117                 pr = &sys->ffparams.iparams[ilist[F_POSRES].iatoms[i]];
3118                 for(d=0; d<DIM; d++)
3119                 {
3120                     if (pr->posres.fcA[d] != 0)
3121                     {
3122                         AbsRef[d] = 1;
3123                     }
3124                 }
3125             }
3126             for(i=0; i<ilist[F_FBPOSRES].nr; i+=2)
3127             {
3128                 /* Check for flat-bottom posres */
3129                 pr = &sys->ffparams.iparams[ilist[F_FBPOSRES].iatoms[i]];
3130                 if (pr->fbposres.k != 0)
3131                 {
3132                     switch(pr->fbposres.geom)
3133                     {
3134                     case efbposresSPHERE:
3135                         AbsRef[XX] = AbsRef[YY] = AbsRef[ZZ] = 1;
3136                         break;
3137                     case efbposresCYLINDER:
3138                         AbsRef[XX] = AbsRef[YY] = 1;
3139                         break;
3140                     case efbposresX: /* d=XX */
3141                     case efbposresY: /* d=YY */
3142                     case efbposresZ: /* d=ZZ */
3143                         d = pr->fbposres.geom - efbposresX;
3144                         AbsRef[d] = 1;
3145                         break;
3146                     default:
3147                         gmx_fatal(FARGS," Invalid geometry for flat-bottom position restraint.\n"
3148                                   "Expected nr between 1 and %d. Found %d\n", efbposresNR-1,
3149                                   pr->fbposres.geom);
3150                     }
3151                 }
3152             }
3153         }
3154     }
3155
3156     return (AbsRef[XX] != 0 && AbsRef[YY] != 0 && AbsRef[ZZ] != 0);
3157 }
3158
3159 void triple_check(const char *mdparin,t_inputrec *ir,gmx_mtop_t *sys,
3160                   warninp_t wi)
3161 {
3162   char err_buf[256];
3163   int  i,m,g,nmol,npct;
3164   gmx_bool bCharge,bAcc;
3165   real gdt_max,*mgrp,mt;
3166   rvec acc;
3167   gmx_mtop_atomloop_block_t aloopb;
3168   gmx_mtop_atomloop_all_t aloop;
3169   t_atom *atom;
3170   ivec AbsRef;
3171   char warn_buf[STRLEN];
3172
3173   set_warning_line(wi,mdparin,-1);
3174
3175   if (EI_DYNAMICS(ir->eI) && !EI_SD(ir->eI) && ir->eI != eiBD &&
3176       ir->comm_mode == ecmNO &&
3177       !(absolute_reference(ir,sys,FALSE,AbsRef) || ir->nsteps <= 10)) {
3178     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");
3179   }
3180
3181     /* Check for pressure coupling with absolute position restraints */
3182     if (ir->epc != epcNO && ir->refcoord_scaling == erscNO)
3183     {
3184         absolute_reference(ir,sys,TRUE,AbsRef);
3185         {
3186             for(m=0; m<DIM; m++)
3187             {
3188                 if (AbsRef[m] && norm2(ir->compress[m]) > 0)
3189                 {
3190                     warning(wi,"You are using pressure coupling with absolute position restraints, this will give artifacts. Use the refcoord_scaling option.");
3191                     break;
3192                 }
3193             }
3194         }
3195     }
3196
3197   bCharge = FALSE;
3198   aloopb = gmx_mtop_atomloop_block_init(sys);
3199   while (gmx_mtop_atomloop_block_next(aloopb,&atom,&nmol)) {
3200     if (atom->q != 0 || atom->qB != 0) {
3201       bCharge = TRUE;
3202     }
3203   }
3204
3205   if (!bCharge) {
3206     if (EEL_FULL(ir->coulombtype)) {
3207       sprintf(err_buf,
3208               "You are using full electrostatics treatment %s for a system without charges.\n"
3209               "This costs a lot of performance for just processing zeros, consider using %s instead.\n",
3210               EELTYPE(ir->coulombtype),EELTYPE(eelCUT));
3211       warning(wi,err_buf);
3212     }
3213   } else {
3214     if (ir->coulombtype == eelCUT && ir->rcoulomb > 0 && !ir->implicit_solvent) {
3215       sprintf(err_buf,
3216               "You are using a plain Coulomb cut-off, which might produce artifacts.\n"
3217               "You might want to consider using %s electrostatics.\n",
3218               EELTYPE(eelPME));
3219       warning_note(wi,err_buf);
3220     }
3221   }
3222
3223   /* Generalized reaction field */
3224   if (ir->opts.ngtc == 0) {
3225     sprintf(err_buf,"No temperature coupling while using coulombtype %s",
3226             eel_names[eelGRF]);
3227     CHECK(ir->coulombtype == eelGRF);
3228   }
3229   else {
3230     sprintf(err_buf,"When using coulombtype = %s"
3231             " ref-t for temperature coupling should be > 0",
3232             eel_names[eelGRF]);
3233     CHECK((ir->coulombtype == eelGRF) && (ir->opts.ref_t[0] <= 0));
3234   }
3235
3236     if (ir->eI == eiSD1 &&
3237         (gmx_mtop_ftype_count(sys,F_CONSTR) > 0 ||
3238          gmx_mtop_ftype_count(sys,F_SETTLE) > 0))
3239     {
3240         sprintf(warn_buf,"With constraints integrator %s is less accurate, consider using %s instead",ei_names[ir->eI],ei_names[eiSD2]);
3241         warning_note(wi,warn_buf);
3242     }
3243
3244   bAcc = FALSE;
3245   for(i=0; (i<sys->groups.grps[egcACC].nr); i++) {
3246     for(m=0; (m<DIM); m++) {
3247       if (fabs(ir->opts.acc[i][m]) > 1e-6) {
3248         bAcc = TRUE;
3249       }
3250     }
3251   }
3252   if (bAcc) {
3253     clear_rvec(acc);
3254     snew(mgrp,sys->groups.grps[egcACC].nr);
3255     aloop = gmx_mtop_atomloop_all_init(sys);
3256     while (gmx_mtop_atomloop_all_next(aloop,&i,&atom)) {
3257       mgrp[ggrpnr(&sys->groups,egcACC,i)] += atom->m;
3258     }
3259     mt = 0.0;
3260     for(i=0; (i<sys->groups.grps[egcACC].nr); i++) {
3261       for(m=0; (m<DIM); m++)
3262         acc[m] += ir->opts.acc[i][m]*mgrp[i];
3263       mt += mgrp[i];
3264     }
3265     for(m=0; (m<DIM); m++) {
3266       if (fabs(acc[m]) > 1e-6) {
3267         const char *dim[DIM] = { "X", "Y", "Z" };
3268         fprintf(stderr,
3269                 "Net Acceleration in %s direction, will %s be corrected\n",
3270                 dim[m],ir->nstcomm != 0 ? "" : "not");
3271         if (ir->nstcomm != 0 && m < ndof_com(ir)) {
3272           acc[m] /= mt;
3273           for (i=0; (i<sys->groups.grps[egcACC].nr); i++)
3274             ir->opts.acc[i][m] -= acc[m];
3275         }
3276       }
3277     }
3278     sfree(mgrp);
3279   }
3280
3281   if (ir->efep != efepNO && ir->fepvals->sc_alpha != 0 &&
3282       !gmx_within_tol(sys->ffparams.reppow,12.0,10*GMX_DOUBLE_EPS)) {
3283     gmx_fatal(FARGS,"Soft-core interactions are only supported with VdW repulsion power 12");
3284   }
3285
3286   if (ir->ePull != epullNO) {
3287     if (ir->pull->grp[0].nat == 0) {
3288         absolute_reference(ir,sys,FALSE,AbsRef);
3289       for(m=0; m<DIM; m++) {
3290         if (ir->pull->dim[m] && !AbsRef[m]) {
3291           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.");
3292           break;
3293         }
3294       }
3295     }
3296
3297     if (ir->pull->eGeom == epullgDIRPBC) {
3298       for(i=0; i<3; i++) {
3299         for(m=0; m<=i; m++) {
3300           if ((ir->epc != epcNO && ir->compress[i][m] != 0) ||
3301               ir->deform[i][m] != 0) {
3302             for(g=1; g<ir->pull->ngrp; g++) {
3303               if (ir->pull->grp[g].vec[m] != 0) {
3304                 gmx_fatal(FARGS,"Can not have dynamic box while using pull geometry '%s' (dim %c)",EPULLGEOM(ir->pull->eGeom),'x'+m);
3305               }
3306             }
3307           }
3308         }
3309       }
3310     }
3311   }
3312
3313   check_disre(sys);
3314 }
3315
3316 void double_check(t_inputrec *ir,matrix box,gmx_bool bConstr,warninp_t wi)
3317 {
3318   real min_size;
3319   gmx_bool bTWIN;
3320   char warn_buf[STRLEN];
3321   const char *ptr;
3322
3323   ptr = check_box(ir->ePBC,box);
3324   if (ptr) {
3325       warning_error(wi,ptr);
3326   }
3327
3328   if (bConstr && ir->eConstrAlg == econtSHAKE) {
3329     if (ir->shake_tol <= 0.0) {
3330       sprintf(warn_buf,"ERROR: shake-tol must be > 0 instead of %g\n",
3331               ir->shake_tol);
3332       warning_error(wi,warn_buf);
3333     }
3334
3335     if (IR_TWINRANGE(*ir) && ir->nstlist > 1) {
3336       sprintf(warn_buf,"With twin-range cut-off's and SHAKE the virial and the pressure are incorrect.");
3337       if (ir->epc == epcNO) {
3338         warning(wi,warn_buf);
3339       } else {
3340           warning_error(wi,warn_buf);
3341       }
3342     }
3343   }
3344
3345   if( (ir->eConstrAlg == econtLINCS) && bConstr) {
3346     /* If we have Lincs constraints: */
3347     if(ir->eI==eiMD && ir->etc==etcNO &&
3348        ir->eConstrAlg==econtLINCS && ir->nLincsIter==1) {
3349       sprintf(warn_buf,"For energy conservation with LINCS, lincs_iter should be 2 or larger.\n");
3350       warning_note(wi,warn_buf);
3351     }
3352
3353     if ((ir->eI == eiCG || ir->eI == eiLBFGS) && (ir->nProjOrder<8)) {
3354       sprintf(warn_buf,"For accurate %s with LINCS constraints, lincs-order should be 8 or more.",ei_names[ir->eI]);
3355       warning_note(wi,warn_buf);
3356     }
3357     if (ir->epc==epcMTTK) {
3358         warning_error(wi,"MTTK not compatible with lincs -- use shake instead.");
3359     }
3360   }
3361
3362   if (ir->LincsWarnAngle > 90.0) {
3363     sprintf(warn_buf,"lincs-warnangle can not be larger than 90 degrees, setting it to 90.\n");
3364     warning(wi,warn_buf);
3365     ir->LincsWarnAngle = 90.0;
3366   }
3367
3368   if (ir->ePBC != epbcNONE) {
3369     if (ir->nstlist == 0) {
3370       warning(wi,"With nstlist=0 atoms are only put into the box at step 0, therefore drifting atoms might cause the simulation to crash.");
3371     }
3372     bTWIN = (ir->rlistlong > ir->rlist);
3373     if (ir->ns_type == ensGRID) {
3374       if (sqr(ir->rlistlong) >= max_cutoff2(ir->ePBC,box)) {
3375           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",
3376                 bTWIN ? (ir->rcoulomb==ir->rlistlong ? "rcoulomb" : "rvdw"):"rlist");
3377           warning_error(wi,warn_buf);
3378       }
3379     } else {
3380       min_size = min(box[XX][XX],min(box[YY][YY],box[ZZ][ZZ]));
3381       if (2*ir->rlistlong >= min_size) {
3382           sprintf(warn_buf,"ERROR: One of the box lengths is smaller than twice the cut-off length. Increase the box size or decrease rlist.");
3383           warning_error(wi,warn_buf);
3384         if (TRICLINIC(box))
3385           fprintf(stderr,"Grid search might allow larger cut-off's than simple search with triclinic boxes.");
3386       }
3387     }
3388   }
3389 }
3390
3391 void check_chargegroup_radii(const gmx_mtop_t *mtop,const t_inputrec *ir,
3392                              rvec *x,
3393                              warninp_t wi)
3394 {
3395     real rvdw1,rvdw2,rcoul1,rcoul2;
3396     char warn_buf[STRLEN];
3397
3398     calc_chargegroup_radii(mtop,x,&rvdw1,&rvdw2,&rcoul1,&rcoul2);
3399
3400     if (rvdw1 > 0)
3401     {
3402         printf("Largest charge group radii for Van der Waals: %5.3f, %5.3f nm\n",
3403                rvdw1,rvdw2);
3404     }
3405     if (rcoul1 > 0)
3406     {
3407         printf("Largest charge group radii for Coulomb:       %5.3f, %5.3f nm\n",
3408                rcoul1,rcoul2);
3409     }
3410
3411     if (ir->rlist > 0)
3412     {
3413         if (rvdw1  + rvdw2  > ir->rlist ||
3414             rcoul1 + rcoul2 > ir->rlist)
3415         {
3416             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);
3417             warning(wi,warn_buf);
3418         }
3419         else
3420         {
3421             /* Here we do not use the zero at cut-off macro,
3422              * since user defined interactions might purposely
3423              * not be zero at the cut-off.
3424              */
3425             if (EVDW_IS_ZERO_AT_CUTOFF(ir->vdwtype) &&
3426                 rvdw1 + rvdw2 > ir->rlist - ir->rvdw)
3427             {
3428                 sprintf(warn_buf,"The sum of the two largest charge group radii (%f) is larger than rlist (%f) - rvdw (%f)\n",
3429                         rvdw1+rvdw2,
3430                         ir->rlist,ir->rvdw);
3431                 if (ir_NVE(ir))
3432                 {
3433                     warning(wi,warn_buf);
3434                 }
3435                 else
3436                 {
3437                     warning_note(wi,warn_buf);
3438                 }
3439             }
3440             if (EEL_IS_ZERO_AT_CUTOFF(ir->coulombtype) &&
3441                 rcoul1 + rcoul2 > ir->rlistlong - ir->rcoulomb)
3442             {
3443                 sprintf(warn_buf,"The sum of the two largest charge group radii (%f) is larger than %s (%f) - rcoulomb (%f)\n",
3444                         rcoul1+rcoul2,
3445                         ir->rlistlong > ir->rlist ? "rlistlong" : "rlist",
3446                         ir->rlistlong,ir->rcoulomb);
3447                 if (ir_NVE(ir))
3448                 {
3449                     warning(wi,warn_buf);
3450                 }
3451                 else
3452                 {
3453                     warning_note(wi,warn_buf);
3454                 }
3455             }
3456         }
3457     }
3458 }