src/kernel/readir.c

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