src/gromacs/gmxpreprocess/readir.c

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