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