src/gromacs/gmxpreprocess/readir.c

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