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