src/kernel/readir.c

   1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
   2  *
   3  *
   4  *                This source code is part of
   5  *
   6  *                 G   R   O   M   A   C   S
   7  *
   8  *          GROningen MAchine for Chemical Simulations
   9  *
  10  *                        VERSION 3.2.0
  11  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  12  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  13  * Copyright (c) 2001-2004, The GROMACS development team,
  14  * check out http://www.gromacs.org for more information.
  15
  16  * This program is free software; you can redistribute it and/or
  17  * modify it under the terms of the GNU General Public License
  18  * as published by the Free Software Foundation; either version 2
  19  * of the License, or (at your option) any later version.
  20  *
  21  * If you want to redistribute modifications, please consider that
  22  * scientific software is very special. Version control is crucial -
  23  * bugs must be traceable. We will be happy to consider code for
  24  * inclusion in the official distribution, but derived work must not
  25  * be called official GROMACS. Details are found in the README & COPYING
  26  * files - if they are missing, get the official version at www.gromacs.org.
  27  *
  28  * To help us fund GROMACS development, we humbly ask that you cite
  29  * the papers on the package - you can find them in the top README file.
  30  *
  31  * For more info, check our website at http://www.gromacs.org
  32  *
  33  * And Hey:
  34  * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
  35  */
  36 #ifdef HAVE_CONFIG_H
  37 #include <config.h>
  38 #endif
  39
  40 #include <ctype.h>
  41 #include <stdlib.h>
  42 #include <limits.h>
  43 #include "sysstuff.h"
  44 #include "smalloc.h"
  45 #include "typedefs.h"
  46 #include "physics.h"
  47 #include "names.h"
  48 #include "gmx_fatal.h"
  49 #include "macros.h"
  50 #include "index.h"
  51 #include "symtab.h"
  52 #include "string2.h"
  53 #include "readinp.h"
  54 #include "warninp.h"
  55 #include "readir.h"
  56 #include "toputil.h"
  57 #include "index.h"
  58 #include "network.h"
  59 #include "vec.h"
  60 #include "pbc.h"
  61 #include "mtop_util.h"
  62 #include "chargegroup.h"
  63 #include "inputrec.h"
  64
  65 #define MAXPTR 254
  66 #define NOGID  255
  67
  68 /* Resource parameters
  69  * Do not change any of these until you read the instruction
  70  * in readinp.h. Some cpp's do not take spaces after the backslash
  71  * (like the c-shell), which will give you a very weird compiler
  72  * message.
  73  */
  74
  75 static char tcgrps[STRLEN],tau_t[STRLEN],ref_t[STRLEN],
  76   acc[STRLEN],accgrps[STRLEN],freeze[STRLEN],frdim[STRLEN],
  77   energy[STRLEN],user1[STRLEN],user2[STRLEN],vcm[STRLEN],xtc_grps[STRLEN],
  78   couple_moltype[STRLEN],orirefitgrp[STRLEN],egptable[STRLEN],egpexcl[STRLEN],
  79   wall_atomtype[STRLEN],wall_density[STRLEN],deform[STRLEN],QMMM[STRLEN];
  80 static char foreign_lambda[STRLEN];
  81 static char **pull_grp;
  82 static char **rot_grp;
  83 static char anneal[STRLEN],anneal_npoints[STRLEN],
  84   anneal_time[STRLEN],anneal_temp[STRLEN];
  85 static char QMmethod[STRLEN],QMbasis[STRLEN],QMcharge[STRLEN],QMmult[STRLEN],
  86   bSH[STRLEN],CASorbitals[STRLEN], CASelectrons[STRLEN],SAon[STRLEN],
  87   SAoff[STRLEN],SAsteps[STRLEN],bTS[STRLEN],bOPT[STRLEN];
  88 static char efield_x[STRLEN],efield_xt[STRLEN],efield_y[STRLEN],
  89   efield_yt[STRLEN],efield_z[STRLEN],efield_zt[STRLEN];
  90
  91 enum {
  92     egrptpALL,         /* All particles have to be a member of a group.     */
  93     egrptpALL_GENREST, /* A rest group with name is generated for particles *
  94                         * that are not part of any group.                   */
  95     egrptpPART,        /* As egrptpALL_GENREST, but no name is generated    *
  96                         * for the rest group.                               */
  97     egrptpONE          /* Merge all selected groups into one group,         *
  98                         * make a rest group for the remaining particles.    */
  99 };
 100
 101
 102 void init_ir(t_inputrec *ir, t_gromppopts *opts)
 103 {
 104   snew(opts->include,STRLEN);
 105   snew(opts->define,STRLEN);
 106 }
 107
 108 static void _low_check(gmx_bool b,char *s,warninp_t wi)
 109 {
 110     if (b)
 111     {
 112         warning_error(wi,s);
 113     }
 114 }
 115
 116 static void check_nst(const char *desc_nst,int nst,
 117                       const char *desc_p,int *p,
 118                       warninp_t wi)
 119 {
 120     char buf[STRLEN];
 121
 122     if (*p > 0 && *p % nst != 0)
 123     {
 124         /* Round up to the next multiple of nst */
 125         *p = ((*p)/nst + 1)*nst;
 126         sprintf(buf,"%s should be a multiple of %s, changing %s to %d\n",
 127                 desc_p,desc_nst,desc_p,*p);
 128         warning(wi,buf);
 129     }
 130 }
 131
 132 static gmx_bool ir_NVE(const t_inputrec *ir)
 133 {
 134     return ((ir->eI == eiMD || EI_VV(ir->eI)) && ir->etc == etcNO);
 135 }
 136
 137 static int lcd(int n1,int n2)
 138 {
 139     int d,i;
 140
 141     d = 1;
 142     for(i=2; (i<=n1 && i<=n2); i++)
 143     {
 144         if (n1 % i == 0 && n2 % i == 0)
 145         {
 146             d = i;
 147         }
 148     }
 149
 150   return d;
 151 }
 152
 153 void check_ir(const char *mdparin,t_inputrec *ir, t_gromppopts *opts,
 154               warninp_t wi)
 155 /* Check internal consistency */
 156 {
 157     /* Strange macro: first one fills the err_buf, and then one can check
 158      * the condition, which will print the message and increase the error
 159      * counter.
 160      */
 161 #define CHECK(b) _low_check(b,err_buf,wi)
 162     char err_buf[256],warn_buf[STRLEN];
 163     int  ns_type=0;
 164     real dt_pcoupl;
 165
 166   set_warning_line(wi,mdparin,-1);
 167
 168   /* BASIC CUT-OFF STUFF */
 169   if (ir->rlist == 0 ||
 170       !((EEL_MIGHT_BE_ZERO_AT_CUTOFF(ir->coulombtype) && ir->rcoulomb > ir->rlist) ||
 171         (EVDW_MIGHT_BE_ZERO_AT_CUTOFF(ir->vdwtype)    && ir->rvdw     > ir->rlist))) {
 172     /* No switched potential and/or no twin-range:
 173      * we can set the long-range cut-off to the maximum of the other cut-offs.
 174      */
 175     ir->rlistlong = max_cutoff(ir->rlist,max_cutoff(ir->rvdw,ir->rcoulomb));
 176   } else if (ir->rlistlong < 0) {
 177     ir->rlistlong = max_cutoff(ir->rlist,max_cutoff(ir->rvdw,ir->rcoulomb));
 178     sprintf(warn_buf,"rlistlong was not set, setting it to %g (no buffer)",
 179             ir->rlistlong);
 180     warning(wi,warn_buf);
 181   }
 182   if (ir->rlistlong == 0 && ir->ePBC != epbcNONE) {
 183       warning_error(wi,"Can not have an infinite cut-off with PBC");
 184   }
 185   if (ir->rlistlong > 0 && (ir->rlist == 0 || ir->rlistlong < ir->rlist)) {
 186       warning_error(wi,"rlistlong can not be shorter than rlist");
 187   }
 188   if (IR_TWINRANGE(*ir) && ir->nstlist <= 0) {
 189       warning_error(wi,"Can not have nstlist<=0 with twin-range interactions");
 190   }
 191
 192     /* GENERAL INTEGRATOR STUFF */
 193     if (!(ir->eI == eiMD || EI_VV(ir->eI)))
 194     {
 195         ir->etc = etcNO;
 196     }
 197     if (!EI_DYNAMICS(ir->eI))
 198     {
 199         ir->epc = epcNO;
 200     }
 201     if (EI_DYNAMICS(ir->eI))
 202     {
 203         if (ir->nstcalcenergy < 0)
 204         {
 205             ir->nstcalcenergy = ir_optimal_nstcalcenergy(ir);
 206             if (ir->nstenergy != 0 && ir->nstenergy < ir->nstcalcenergy)
 207             {
 208                 /* nstcalcenergy larger than nstener does not make sense.
 209                  * We ideally want nstcalcenergy=nstener.
 210                  */
 211                 if (ir->nstlist > 0)
 212                 {
 213                     ir->nstcalcenergy = lcd(ir->nstenergy,ir->nstlist);
 214                 }
 215                 else
 216                 {
 217                     ir->nstcalcenergy = ir->nstenergy;
 218                 }
 219             }
 220         }
 221         if (ir->epc != epcNO)
 222         {
 223             if (ir->nstpcouple < 0)
 224             {
 225                 ir->nstpcouple = ir_optimal_nstpcouple(ir);
 226             }
 227         }
 228         if (IR_TWINRANGE(*ir))
 229         {
 230             check_nst("nstlist",ir->nstlist,
 231                       "nstcalcenergy",&ir->nstcalcenergy,wi);
 232             if (ir->epc != epcNO)
 233             {
 234                 check_nst("nstlist",ir->nstlist,
 235                           "nstpcouple",&ir->nstpcouple,wi);
 236             }
 237         }
 238
 239         if (ir->nstcalcenergy > 1)
 240         {
 241             /* for storing exact averages nstenergy should be
 242              * a multiple of nstcalcenergy
 243              */
 244             check_nst("nstcalcenergy",ir->nstcalcenergy,
 245                       "nstenergy",&ir->nstenergy,wi);
 246             if (ir->efep != efepNO)
 247             {
 248                 /* nstdhdl should be a multiple of nstcalcenergy */
 249                 check_nst("nstcalcenergy",ir->nstcalcenergy,
 250                           "nstdhdl",&ir->nstdhdl,wi);
 251             }
 252         }
 253     }
 254
 255   /* LD STUFF */
 256   if ((EI_SD(ir->eI) || ir->eI == eiBD) &&
 257       ir->bContinuation && ir->ld_seed != -1) {
 258       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)");
 259   }
 260
 261   /* TPI STUFF */
 262   if (EI_TPI(ir->eI)) {
 263     sprintf(err_buf,"TPI only works with pbc = %s",epbc_names[epbcXYZ]);
 264     CHECK(ir->ePBC != epbcXYZ);
 265     sprintf(err_buf,"TPI only works with ns = %s",ens_names[ensGRID]);
 266     CHECK(ir->ns_type != ensGRID);
 267     sprintf(err_buf,"with TPI nstlist should be larger than zero");
 268     CHECK(ir->nstlist <= 0);
 269     sprintf(err_buf,"TPI does not work with full electrostatics other than PME");
 270     CHECK(EEL_FULL(ir->coulombtype) && !EEL_PME(ir->coulombtype));
 271   }
 272
 273   /* SHAKE / LINCS */
 274   if ( (opts->nshake > 0) && (opts->bMorse) ) {
 275     sprintf(warn_buf,
 276             "Using morse bond-potentials while constraining bonds is useless");
 277     warning(wi,warn_buf);
 278   }
 279
 280   sprintf(err_buf,"shake_tol must be > 0 instead of %g while using shake",
 281           ir->shake_tol);
 282   CHECK(((ir->shake_tol <= 0.0) && (opts->nshake>0) &&
 283          (ir->eConstrAlg == econtSHAKE)));
 284
 285   /* PBC/WALLS */
 286   sprintf(err_buf,"walls only work with pbc=%s",epbc_names[epbcXY]);
 287   CHECK(ir->nwall && ir->ePBC!=epbcXY);
 288
 289   /* VACUUM STUFF */
 290   if (ir->ePBC != epbcXYZ && ir->nwall != 2) {
 291     if (ir->ePBC == epbcNONE) {
 292       if (ir->epc != epcNO) {
 293           warning(wi,"Turning off pressure coupling for vacuum system");
 294           ir->epc = epcNO;
 295       }
 296     } else {
 297       sprintf(err_buf,"Can not have pressure coupling with pbc=%s",
 298               epbc_names[ir->ePBC]);
 299       CHECK(ir->epc != epcNO);
 300     }
 301     sprintf(err_buf,"Can not have Ewald with pbc=%s",epbc_names[ir->ePBC]);
 302     CHECK(EEL_FULL(ir->coulombtype));
 303
 304     sprintf(err_buf,"Can not have dispersion correction with pbc=%s",
 305             epbc_names[ir->ePBC]);
 306     CHECK(ir->eDispCorr != edispcNO);
 307   }
 308
 309   if (ir->rlist == 0.0) {
 310     sprintf(err_buf,"can only have neighborlist cut-off zero (=infinite)\n"
 311             "with coulombtype = %s or coulombtype = %s\n"
 312             "without periodic boundary conditions (pbc = %s) and\n"
 313             "rcoulomb and rvdw set to zero",
 314             eel_names[eelCUT],eel_names[eelUSER],epbc_names[epbcNONE]);
 315     CHECK(((ir->coulombtype != eelCUT) && (ir->coulombtype != eelUSER)) ||
 316           (ir->ePBC     != epbcNONE) ||
 317           (ir->rcoulomb != 0.0)      || (ir->rvdw != 0.0));
 318
 319     if (ir->nstlist < 0) {
 320         warning_error(wi,"Can not have heuristic neighborlist updates without cut-off");
 321     }
 322     if (ir->nstlist > 0) {
 323         warning_note(wi,"Simulating without cut-offs is usually (slightly) faster with nstlist=0, nstype=simple and particle decomposition");
 324     }
 325   }
 326
 327   /* COMM STUFF */
 328   if (ir->nstcomm == 0) {
 329     ir->comm_mode = ecmNO;
 330   }
 331   if (ir->comm_mode != ecmNO) {
 332     if (ir->nstcomm < 0) {
 333         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");
 334       ir->nstcomm = abs(ir->nstcomm);
 335     }
 336
 337     if (ir->nstcalcenergy > 0 && ir->nstcomm < ir->nstcalcenergy) {
 338         warning_note(wi,"nstcomm < nstcalcenergy defeats the purpose of nstcalcenergy, setting nstcomm to nstcalcenergy");
 339       ir->nstcomm = ir->nstcalcenergy;
 340     }
 341
 342     if (ir->comm_mode == ecmANGULAR) {
 343       sprintf(err_buf,"Can not remove the rotation around the center of mass with periodic molecules");
 344       CHECK(ir->bPeriodicMols);
 345       if (ir->ePBC != epbcNONE)
 346           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).");
 347     }
 348   }
 349
 350   if (EI_STATE_VELOCITY(ir->eI) && ir->ePBC == epbcNONE && ir->comm_mode != ecmANGULAR) {
 351       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.");
 352   }
 353
 354   sprintf(err_buf,"Free-energy not implemented for Ewald and PPPM");
 355   CHECK((ir->coulombtype==eelEWALD || ir->coulombtype==eelPPPM)
 356         && (ir->efep!=efepNO));
 357
 358   sprintf(err_buf,"Twin-range neighbour searching (NS) with simple NS"
 359           " algorithm not implemented");
 360   CHECK(((ir->rcoulomb > ir->rlist) || (ir->rvdw > ir->rlist))
 361         && (ir->ns_type == ensSIMPLE));
 362
 363     /* TEMPERATURE COUPLING */
 364     if (ir->etc == etcYES)
 365     {
 366         ir->etc = etcBERENDSEN;
 367         warning_note(wi,"Old option for temperature coupling given: "
 368                      "changing \"yes\" to \"Berendsen\"\n");
 369     }
 370
 371     if (ir->etc == etcNOSEHOOVER)
 372     {
 373         if (ir->opts.nhchainlength < 1)
 374         {
 375             sprintf(warn_buf,"number of Nose-Hoover chains (currently %d) cannot be less than 1,reset to 1\n",ir->opts.nhchainlength);
 376             ir->opts.nhchainlength =1;
 377             warning(wi,warn_buf);
 378         }
 379
 380         if (ir->etc==etcNOSEHOOVER && !EI_VV(ir->eI) && ir->opts.nhchainlength > 1)
 381         {
 382             warning_note(wi,"leapfrog does not yet support Nose-Hoover chains, nhchainlength reset to 1");
 383             ir->opts.nhchainlength = 1;
 384         }
 385     }
 386     else
 387     {
 388         ir->opts.nhchainlength = 0;
 389     }
 390
 391     if (ir->etc == etcBERENDSEN)
 392     {
 393         sprintf(warn_buf,"The %s thermostat does not generate the correct kinetic energy distribution. You might want to consider using the %s thermostat.",
 394                 ETCOUPLTYPE(ir->etc),ETCOUPLTYPE(etcVRESCALE));
 395         warning_note(wi,warn_buf);
 396     }
 397
 398     if ((ir->etc==etcNOSEHOOVER || ir->etc==etcANDERSEN || ir->etc==etcANDERSENINTERVAL)
 399         && ir->epc==epcBERENDSEN)
 400     {
 401         sprintf(warn_buf,"Using Berendsen pressure coupling invalidates the "
 402                 "true ensemble for the thermostat");
 403         warning(wi,warn_buf);
 404     }
 405
 406     /* PRESSURE COUPLING */
 407     if (ir->epc == epcISOTROPIC)
 408     {
 409         ir->epc = epcBERENDSEN;
 410         warning_note(wi,"Old option for pressure coupling given: "
 411                      "changing \"Isotropic\" to \"Berendsen\"\n");
 412     }
 413
 414     if (ir->epc != epcNO)
 415     {
 416         dt_pcoupl = ir->nstpcouple*ir->delta_t;
 417
 418         sprintf(err_buf,"tau-p must be > 0 instead of %g\n",ir->tau_p);
 419         CHECK(ir->tau_p <= 0);
 420
 421         if (ir->tau_p/dt_pcoupl < pcouple_min_integration_steps(ir->epc))
 422         {
 423             sprintf(warn_buf,"For proper integration of the %s barostat, tau-p (%g) should be at least %d times larger than nstpcouple*dt (%g)",
 424                     EPCOUPLTYPE(ir->epc),ir->tau_p,pcouple_min_integration_steps(ir->epc),dt_pcoupl);
 425             warning(wi,warn_buf);
 426         }
 427
 428         sprintf(err_buf,"compressibility must be > 0 when using pressure"
 429                 " coupling %s\n",EPCOUPLTYPE(ir->epc));
 430         CHECK(ir->compress[XX][XX] < 0 || ir->compress[YY][YY] < 0 ||
 431               ir->compress[ZZ][ZZ] < 0 ||
 432               (trace(ir->compress) == 0 && ir->compress[YY][XX] <= 0 &&
 433                ir->compress[ZZ][XX] <= 0 && ir->compress[ZZ][YY] <= 0));
 434
 435         sprintf(err_buf,"pressure coupling with PPPM not implemented, use PME");
 436         CHECK(ir->coulombtype == eelPPPM);
 437
 438     }
 439     else if (ir->coulombtype == eelPPPM)
 440     {
 441         sprintf(warn_buf,"The pressure with PPPM is incorrect, if you need the pressure use PME");
 442         warning(wi,warn_buf);
 443     }
 444
 445     if (EI_VV(ir->eI))
 446     {
 447         if (ir->epc > epcNO)
 448         {
 449             if (ir->epc!=epcMTTK)
 450             {
 451                 warning_error(wi,"NPT only defined for vv using Martyna-Tuckerman-Tobias-Klein equations");
 452             }
 453         }
 454     }
 455
 456   /* ELECTROSTATICS */
 457   /* More checks are in triple check (grompp.c) */
 458     if (ir->coulombtype == eelPPPM)
 459     {
 460         warning_error(wi,"PPPM is not functional in the current version, we plan to implement PPPM through a small modification of the PME code");
 461     }
 462
 463   if (ir->coulombtype == eelSWITCH) {
 464     sprintf(warn_buf,"coulombtype = %s is only for testing purposes and can lead to serious artifacts, advice: use coulombtype = %s",
 465             eel_names[ir->coulombtype],
 466             eel_names[eelRF_ZERO]);
 467     warning(wi,warn_buf);
 468   }
 469
 470   if (ir->epsilon_r!=1 && ir->implicit_solvent==eisGBSA) {
 471     sprintf(warn_buf,"epsilon-r = %g with GB implicit solvent, will use this value for inner dielectric",ir->epsilon_r);
 472     warning_note(wi,warn_buf);
 473   }
 474
 475   if (EEL_RF(ir->coulombtype) && ir->epsilon_rf==1 && ir->epsilon_r!=1) {
 476     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);
 477     warning(wi,warn_buf);
 478     ir->epsilon_rf = ir->epsilon_r;
 479     ir->epsilon_r  = 1.0;
 480   }
 481
 482   if (getenv("GALACTIC_DYNAMICS") == NULL) {
 483     sprintf(err_buf,"epsilon-r must be >= 0 instead of %g\n",ir->epsilon_r);
 484     CHECK(ir->epsilon_r < 0);
 485   }
 486
 487   if (EEL_RF(ir->coulombtype)) {
 488     /* reaction field (at the cut-off) */
 489
 490     if (ir->coulombtype == eelRF_ZERO) {
 491        sprintf(err_buf,"With coulombtype = %s, epsilon-rf must be 0",
 492                eel_names[ir->coulombtype]);
 493       CHECK(ir->epsilon_rf != 0);
 494     }
 495
 496     sprintf(err_buf,"epsilon-rf must be >= epsilon-r");
 497     CHECK((ir->epsilon_rf < ir->epsilon_r && ir->epsilon_rf != 0) ||
 498           (ir->epsilon_r == 0));
 499     if (ir->epsilon_rf == ir->epsilon_r) {
 500       sprintf(warn_buf,"Using epsilon-rf = epsilon-r with %s does not make sense",
 501               eel_names[ir->coulombtype]);
 502       warning(wi,warn_buf);
 503     }
 504   }
 505   /* Allow rlist>rcoulomb for tabulated long range stuff. This just
 506    * means the interaction is zero outside rcoulomb, but it helps to
 507    * provide accurate energy conservation.
 508    */
 509   if (EEL_MIGHT_BE_ZERO_AT_CUTOFF(ir->coulombtype)) {
 510     if (EEL_SWITCHED(ir->coulombtype)) {
 511       sprintf(err_buf,
 512               "With coulombtype = %s rcoulomb_switch must be < rcoulomb",
 513               eel_names[ir->coulombtype]);
 514       CHECK(ir->rcoulomb_switch >= ir->rcoulomb);
 515     }
 516   } else if (ir->coulombtype == eelCUT || EEL_RF(ir->coulombtype)) {
 517     sprintf(err_buf,"With coulombtype = %s, rcoulomb must be >= rlist",
 518             eel_names[ir->coulombtype]);
 519     CHECK(ir->rlist > ir->rcoulomb);
 520   }
 521
 522   if (EEL_FULL(ir->coulombtype)) {
 523     if (ir->coulombtype==eelPMESWITCH || ir->coulombtype==eelPMEUSER ||
 524         ir->coulombtype==eelPMEUSERSWITCH) {
 525       sprintf(err_buf,"With coulombtype = %s, rcoulomb must be <= rlist",
 526               eel_names[ir->coulombtype]);
 527       CHECK(ir->rcoulomb > ir->rlist);
 528     } else {
 529       if (ir->coulombtype == eelPME) {
 530         sprintf(err_buf,
 531                 "With coulombtype = %s, rcoulomb must be equal to rlist\n"
 532                 "If you want optimal energy conservation or exact integration use %s",
 533                 eel_names[ir->coulombtype],eel_names[eelPMESWITCH]);
 534       } else {
 535         sprintf(err_buf,
 536                 "With coulombtype = %s, rcoulomb must be equal to rlist",
 537                 eel_names[ir->coulombtype]);
 538       }
 539       CHECK(ir->rcoulomb != ir->rlist);
 540     }
 541   }
 542
 543   if (EEL_PME(ir->coulombtype)) {
 544     if (ir->pme_order < 3) {
 545         warning_error(wi,"pme-order can not be smaller than 3");
 546     }
 547   }
 548
 549   if (ir->nwall==2 && EEL_FULL(ir->coulombtype)) {
 550     if (ir->ewald_geometry == eewg3D) {
 551       sprintf(warn_buf,"With pbc=%s you should use ewald-geometry=%s",
 552               epbc_names[ir->ePBC],eewg_names[eewg3DC]);
 553       warning(wi,warn_buf);
 554     }
 555     /* This check avoids extra pbc coding for exclusion corrections */
 556     sprintf(err_buf,"wall-ewald-zfac should be >= 2");
 557     CHECK(ir->wall_ewald_zfac < 2);
 558   }
 559
 560   if (EVDW_SWITCHED(ir->vdwtype)) {
 561     sprintf(err_buf,"With vdwtype = %s rvdw-switch must be < rvdw",
 562             evdw_names[ir->vdwtype]);
 563     CHECK(ir->rvdw_switch >= ir->rvdw);
 564   } else if (ir->vdwtype == evdwCUT) {
 565     sprintf(err_buf,"With vdwtype = %s, rvdw must be >= rlist",evdw_names[ir->vdwtype]);
 566     CHECK(ir->rlist > ir->rvdw);
 567   }
 568   if (EEL_IS_ZERO_AT_CUTOFF(ir->coulombtype)
 569       && (ir->rlistlong <= ir->rcoulomb)) {
 570     sprintf(warn_buf,"For energy conservation with switch/shift potentials, %s should be 0.1 to 0.3 nm larger than rcoulomb.",
 571             IR_TWINRANGE(*ir) ? "rlistlong" : "rlist");
 572     warning_note(wi,warn_buf);
 573   }
 574   if (EVDW_SWITCHED(ir->vdwtype) && (ir->rlistlong <= ir->rvdw)) {
 575     sprintf(warn_buf,"For energy conservation with switch/shift potentials, %s should be 0.1 to 0.3 nm larger than rvdw.",
 576             IR_TWINRANGE(*ir) ? "rlistlong" : "rlist");
 577     warning_note(wi,warn_buf);
 578   }
 579
 580   if (ir->vdwtype == evdwUSER && ir->eDispCorr != edispcNO) {
 581       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.");
 582   }
 583
 584   if (ir->nstlist == -1) {
 585     sprintf(err_buf,
 586             "nstlist=-1 only works with switched or shifted potentials,\n"
 587             "suggestion: use vdw-type=%s and coulomb-type=%s",
 588             evdw_names[evdwSHIFT],eel_names[eelPMESWITCH]);
 589     CHECK(!(EEL_MIGHT_BE_ZERO_AT_CUTOFF(ir->coulombtype) &&
 590             EVDW_MIGHT_BE_ZERO_AT_CUTOFF(ir->vdwtype)));
 591
 592     sprintf(err_buf,"With nstlist=-1 rvdw and rcoulomb should be smaller than rlist to account for diffusion and possibly charge-group radii");
 593     CHECK(ir->rvdw >= ir->rlist || ir->rcoulomb >= ir->rlist);
 594   }
 595   sprintf(err_buf,"nstlist can not be smaller than -1");
 596   CHECK(ir->nstlist < -1);
 597
 598   if (ir->eI == eiLBFGS && (ir->coulombtype==eelCUT || ir->vdwtype==evdwCUT)
 599      && ir->rvdw != 0) {
 600     warning(wi,"For efficient BFGS minimization, use switch/shift/pme instead of cut-off.");
 601   }
 602
 603   if (ir->eI == eiLBFGS && ir->nbfgscorr <= 0) {
 604     warning(wi,"Using L-BFGS with nbfgscorr<=0 just gets you steepest descent.");
 605   }
 606
 607   /* FREE ENERGY */
 608   if (ir->efep != efepNO) {
 609     sprintf(err_buf,"The soft-core power is %d and can only be 1 or 2",
 610             ir->sc_power);
 611     CHECK(ir->sc_alpha!=0 && ir->sc_power!=1 && ir->sc_power!=2);
 612   }
 613
 614     /* ENERGY CONSERVATION */
 615     if (ir_NVE(ir))
 616     {
 617         if (!EVDW_MIGHT_BE_ZERO_AT_CUTOFF(ir->vdwtype) && ir->rvdw > 0)
 618         {
 619             sprintf(warn_buf,"You are using a cut-off for VdW interactions with NVE, for good energy conservation use vdwtype = %s (possibly with DispCorr)",
 620                     evdw_names[evdwSHIFT]);
 621             warning_note(wi,warn_buf);
 622         }
 623         if (!EEL_MIGHT_BE_ZERO_AT_CUTOFF(ir->coulombtype) && ir->rcoulomb > 0)
 624         {
 625             sprintf(warn_buf,"You are using a cut-off for electrostatics with NVE, for good energy conservation use coulombtype = %s or %s",
 626                     eel_names[eelPMESWITCH],eel_names[eelRF_ZERO]);
 627             warning_note(wi,warn_buf);
 628         }
 629     }
 630
 631   /* IMPLICIT SOLVENT */
 632   if(ir->coulombtype==eelGB_NOTUSED)
 633   {
 634     ir->coulombtype=eelCUT;
 635     ir->implicit_solvent=eisGBSA;
 636     fprintf(stderr,"Note: Old option for generalized born electrostatics given:\n"
 637             "Changing coulombtype from \"generalized-born\" to \"cut-off\" and instead\n"
 638             "setting implicit-solvent value to \"GBSA\" in input section.\n");
 639   }
 640
 641   if(ir->sa_algorithm==esaSTILL)
 642   {
 643     sprintf(err_buf,"Still SA algorithm not available yet, use %s or %s instead\n",esa_names[esaAPPROX],esa_names[esaNO]);
 644     CHECK(ir->sa_algorithm == esaSTILL);
 645   }
 646
 647   if(ir->implicit_solvent==eisGBSA)
 648   {
 649     sprintf(err_buf,"With GBSA implicit solvent, rgbradii must be equal to rlist.");
 650     CHECK(ir->rgbradii != ir->rlist);
 651
 652     if(ir->coulombtype!=eelCUT)
 653           {
 654                   sprintf(err_buf,"With GBSA, coulombtype must be equal to %s\n",eel_names[eelCUT]);
 655                   CHECK(ir->coulombtype!=eelCUT);
 656           }
 657           if(ir->vdwtype!=evdwCUT)
 658           {
 659                   sprintf(err_buf,"With GBSA, vdw-type must be equal to %s\n",evdw_names[evdwCUT]);
 660                   CHECK(ir->vdwtype!=evdwCUT);
 661           }
 662     if(ir->nstgbradii<1)
 663     {
 664       sprintf(warn_buf,"Using GBSA with nstgbradii<1, setting nstgbradii=1");
 665       warning_note(wi,warn_buf);
 666       ir->nstgbradii=1;
 667     }
 668     if(ir->sa_algorithm==esaNO)
 669     {
 670       sprintf(warn_buf,"No SA (non-polar) calculation requested together with GB. Are you sure this is what you want?\n");
 671       warning_note(wi,warn_buf);
 672     }
 673     if(ir->sa_surface_tension<0 && ir->sa_algorithm!=esaNO)
 674     {
 675       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");
 676       warning_note(wi,warn_buf);
 677
 678       if(ir->gb_algorithm==egbSTILL)
 679       {
 680         ir->sa_surface_tension = 0.0049 * CAL2JOULE * 100;
 681       }
 682       else
 683       {
 684         ir->sa_surface_tension = 0.0054 * CAL2JOULE * 100;
 685       }
 686     }
 687     if(ir->sa_surface_tension==0 && ir->sa_algorithm!=esaNO)
 688     {
 689       sprintf(err_buf, "Surface tension set to 0 while SA-calculation requested\n");
 690       CHECK(ir->sa_surface_tension==0 && ir->sa_algorithm!=esaNO);
 691     }
 692
 693   }
 694 }
 695
 696 static int str_nelem(const char *str,int maxptr,char *ptr[])
 697 {
 698   int  np=0;
 699   char *copy0,*copy;
 700
 701   copy0=strdup(str);
 702   copy=copy0;
 703   ltrim(copy);
 704   while (*copy != '\0') {
 705     if (np >= maxptr)
 706       gmx_fatal(FARGS,"Too many groups on line: '%s' (max is %d)",
 707                   str,maxptr);
 708     if (ptr)
 709       ptr[np]=copy;
 710     np++;
 711     while ((*copy != '\0') && !isspace(*copy))
 712       copy++;
 713     if (*copy != '\0') {
 714       *copy='\0';
 715       copy++;
 716     }
 717     ltrim(copy);
 718   }
 719   if (ptr == NULL)
 720     sfree(copy0);
 721
 722   return np;
 723 }
 724
 725 static void parse_n_double(char *str,int *n,double **r)
 726 {
 727   char *ptr[MAXPTR];
 728   int  i;
 729
 730   *n = str_nelem(str,MAXPTR,ptr);
 731
 732   snew(*r,*n);
 733   for(i=0; i<*n; i++) {
 734     (*r)[i] = strtod(ptr[i],NULL);
 735   }
 736 }
 737
 738 static void do_wall_params(t_inputrec *ir,
 739                            char *wall_atomtype, char *wall_density,
 740                            t_gromppopts *opts)
 741 {
 742     int  nstr,i;
 743     char *names[MAXPTR];
 744     double dbl;
 745
 746     opts->wall_atomtype[0] = NULL;
 747     opts->wall_atomtype[1] = NULL;
 748
 749     ir->wall_atomtype[0] = -1;
 750     ir->wall_atomtype[1] = -1;
 751     ir->wall_density[0] = 0;
 752     ir->wall_density[1] = 0;
 753
 754     if (ir->nwall > 0)
 755     {
 756         nstr = str_nelem(wall_atomtype,MAXPTR,names);
 757         if (nstr != ir->nwall)
 758         {
 759             gmx_fatal(FARGS,"Expected %d elements for wall_atomtype, found %d",
 760                       ir->nwall,nstr);
 761         }
 762         for(i=0; i<ir->nwall; i++)
 763         {
 764             opts->wall_atomtype[i] = strdup(names[i]);
 765         }
 766
 767         if (ir->wall_type == ewt93 || ir->wall_type == ewt104) {
 768             nstr = str_nelem(wall_density,MAXPTR,names);
 769             if (nstr != ir->nwall)
 770             {
 771                 gmx_fatal(FARGS,"Expected %d elements for wall-density, found %d",ir->nwall,nstr);
 772             }
 773             for(i=0; i<ir->nwall; i++)
 774             {
 775                 sscanf(names[i],"%lf",&dbl);
 776                 if (dbl <= 0)
 777                 {
 778                     gmx_fatal(FARGS,"wall-density[%d] = %f\n",i,dbl);
 779                 }
 780                 ir->wall_density[i] = dbl;
 781             }
 782         }
 783     }
 784 }
 785
 786 static void add_wall_energrps(gmx_groups_t *groups,int nwall,t_symtab *symtab)
 787 {
 788   int  i;
 789   t_grps *grps;
 790   char str[STRLEN];
 791
 792   if (nwall > 0) {
 793     srenew(groups->grpname,groups->ngrpname+nwall);
 794     grps = &(groups->grps[egcENER]);
 795     srenew(grps->nm_ind,grps->nr+nwall);
 796     for(i=0; i<nwall; i++) {
 797       sprintf(str,"wall%d",i);
 798       groups->grpname[groups->ngrpname] = put_symtab(symtab,str);
 799       grps->nm_ind[grps->nr++] = groups->ngrpname++;
 800     }
 801   }
 802 }
 803
 804 void get_ir(const char *mdparin,const char *mdparout,
 805             t_inputrec *ir,t_gromppopts *opts,
 806             warninp_t wi)
 807 {
 808   char      *dumstr[2];
 809   double    dumdub[2][6];
 810   t_inpfile *inp;
 811   const char *tmp;
 812   int       i,j,m,ninp;
 813   char      warn_buf[STRLEN];
 814
 815   inp = read_inpfile(mdparin, &ninp, NULL, wi);
 816
 817   snew(dumstr[0],STRLEN);
 818   snew(dumstr[1],STRLEN);
 819
 820   REM_TYPE("title");
 821   REM_TYPE("cpp");
 822   REM_TYPE("domain-decomposition");
 823   REPL_TYPE("unconstrained-start","continuation");
 824   REM_TYPE("dihre-tau");
 825   REM_TYPE("nstdihreout");
 826   REM_TYPE("nstcheckpoint");
 827
 828   CCTYPE ("VARIOUS PREPROCESSING OPTIONS");
 829   CTYPE ("Preprocessor information: use cpp syntax.");
 830   CTYPE ("e.g.: -I/home/joe/doe -I/home/mary/roe");
 831   STYPE ("include",     opts->include,  NULL);
 832   CTYPE ("e.g.: -DPOSRES -DFLEXIBLE (note these variable names are case sensitive)");
 833   STYPE ("define",      opts->define,   NULL);
 834
 835   CCTYPE ("RUN CONTROL PARAMETERS");
 836   EETYPE("integrator",  ir->eI,         ei_names);
 837   CTYPE ("Start time and timestep in ps");
 838   RTYPE ("tinit",       ir->init_t,     0.0);
 839   RTYPE ("dt",          ir->delta_t,    0.001);
 840   STEPTYPE ("nsteps",   ir->nsteps,     0);
 841   CTYPE ("For exact run continuation or redoing part of a run");
 842   STEPTYPE ("init-step",ir->init_step,  0);
 843   CTYPE ("Part index is updated automatically on checkpointing (keeps files separate)");
 844   ITYPE ("simulation-part", ir->simulation_part, 1);
 845   CTYPE ("mode for center of mass motion removal");
 846   EETYPE("comm-mode",   ir->comm_mode,  ecm_names);
 847   CTYPE ("number of steps for center of mass motion removal");
 848   ITYPE ("nstcomm",     ir->nstcomm,    10);
 849   CTYPE ("group(s) for center of mass motion removal");
 850   STYPE ("comm-grps",   vcm,            NULL);
 851
 852   CCTYPE ("LANGEVIN DYNAMICS OPTIONS");
 853   CTYPE ("Friction coefficient (amu/ps) and random seed");
 854   RTYPE ("bd-fric",     ir->bd_fric,    0.0);
 855   ITYPE ("ld-seed",     ir->ld_seed,    1993);
 856
 857   /* Em stuff */
 858   CCTYPE ("ENERGY MINIMIZATION OPTIONS");
 859   CTYPE ("Force tolerance and initial step-size");
 860   RTYPE ("emtol",       ir->em_tol,     10.0);
 861   RTYPE ("emstep",      ir->em_stepsize,0.01);
 862   CTYPE ("Max number of iterations in relax-shells");
 863   ITYPE ("niter",       ir->niter,      20);
 864   CTYPE ("Step size (ps^2) for minimization of flexible constraints");
 865   RTYPE ("fcstep",      ir->fc_stepsize, 0);
 866   CTYPE ("Frequency of steepest descents steps when doing CG");
 867   ITYPE ("nstcgsteep",  ir->nstcgsteep, 1000);
 868   ITYPE ("nbfgscorr",   ir->nbfgscorr,  10);
 869
 870   CCTYPE ("TEST PARTICLE INSERTION OPTIONS");
 871   RTYPE ("rtpi",        ir->rtpi,       0.05);
 872
 873   /* Output options */
 874   CCTYPE ("OUTPUT CONTROL OPTIONS");
 875   CTYPE ("Output frequency for coords (x), velocities (v) and forces (f)");
 876   ITYPE ("nstxout",     ir->nstxout,    100);
 877   ITYPE ("nstvout",     ir->nstvout,    100);
 878   ITYPE ("nstfout",     ir->nstfout,    0);
 879   ir->nstcheckpoint = 1000;
 880   CTYPE ("Output frequency for energies to log file and energy file");
 881   ITYPE ("nstlog",      ir->nstlog,     100);
 882   ITYPE ("nstcalcenergy",ir->nstcalcenergy,     -1);
 883   ITYPE ("nstenergy",   ir->nstenergy,  100);
 884   CTYPE ("Output frequency and precision for .xtc file");
 885   ITYPE ("nstxtcout",   ir->nstxtcout,  0);
 886   RTYPE ("xtc-precision",ir->xtcprec,   1000.0);
 887   CTYPE ("This selects the subset of atoms for the .xtc file. You can");
 888   CTYPE ("select multiple groups. By default all atoms will be written.");
 889   STYPE ("xtc-grps",    xtc_grps,       NULL);
 890   CTYPE ("Selection of energy groups");
 891   STYPE ("energygrps",  energy,         NULL);
 892
 893   /* Neighbor searching */
 894   CCTYPE ("NEIGHBORSEARCHING PARAMETERS");
 895   CTYPE ("nblist update frequency");
 896   ITYPE ("nstlist",     ir->nstlist,    10);
 897   CTYPE ("ns algorithm (simple or grid)");
 898   EETYPE("ns-type",     ir->ns_type,    ens_names);
 899   /* set ndelta to the optimal value of 2 */
 900   ir->ndelta = 2;
 901   CTYPE ("Periodic boundary conditions: xyz, no, xy");
 902   EETYPE("pbc",         ir->ePBC,       epbc_names);
 903   EETYPE("periodic-molecules", ir->bPeriodicMols, yesno_names);
 904   CTYPE ("nblist cut-off");
 905   RTYPE ("rlist",       ir->rlist,      1.0);
 906   CTYPE ("long-range cut-off for switched potentials");
 907   RTYPE ("rlistlong",   ir->rlistlong,  -1);
 908
 909   /* Electrostatics */
 910   CCTYPE ("OPTIONS FOR ELECTROSTATICS AND VDW");
 911   CTYPE ("Method for doing electrostatics");
 912   EETYPE("coulombtype", ir->coulombtype,    eel_names);
 913   CTYPE ("cut-off lengths");
 914   RTYPE ("rcoulomb-switch",     ir->rcoulomb_switch,    0.0);
 915   RTYPE ("rcoulomb",    ir->rcoulomb,   1.0);
 916   CTYPE ("Relative dielectric constant for the medium and the reaction field");
 917   RTYPE ("epsilon-r",   ir->epsilon_r,  1.0);
 918   RTYPE ("epsilon-rf",  ir->epsilon_rf, 1.0);
 919   CTYPE ("Method for doing Van der Waals");
 920   EETYPE("vdw-type",    ir->vdwtype,    evdw_names);
 921   CTYPE ("cut-off lengths");
 922   RTYPE ("rvdw-switch", ir->rvdw_switch,        0.0);
 923   RTYPE ("rvdw",        ir->rvdw,       1.0);
 924   CTYPE ("Apply long range dispersion corrections for Energy and Pressure");
 925   EETYPE("DispCorr",    ir->eDispCorr,  edispc_names);
 926   CTYPE ("Extension of the potential lookup tables beyond the cut-off");
 927   RTYPE ("table-extension", ir->tabext, 1.0);
 928   CTYPE ("Seperate tables between energy group pairs");
 929   STYPE ("energygrp-table", egptable,   NULL);
 930   CTYPE ("Spacing for the PME/PPPM FFT grid");
 931   RTYPE ("fourierspacing", opts->fourierspacing,0.12);
 932   CTYPE ("FFT grid size, when a value is 0 fourierspacing will be used");
 933   ITYPE ("fourier-nx",  ir->nkx,         0);
 934   ITYPE ("fourier-ny",  ir->nky,         0);
 935   ITYPE ("fourier-nz",  ir->nkz,         0);
 936   CTYPE ("EWALD/PME/PPPM parameters");
 937   ITYPE ("pme-order",   ir->pme_order,   4);
 938   RTYPE ("ewald-rtol",  ir->ewald_rtol, 0.00001);
 939   EETYPE("ewald-geometry", ir->ewald_geometry, eewg_names);
 940   RTYPE ("epsilon-surface", ir->epsilon_surface, 0.0);
 941   EETYPE("optimize-fft",ir->bOptFFT,  yesno_names);
 942
 943   CCTYPE("IMPLICIT SOLVENT ALGORITHM");
 944   EETYPE("implicit-solvent", ir->implicit_solvent, eis_names);
 945
 946   CCTYPE ("GENERALIZED BORN ELECTROSTATICS");
 947   CTYPE ("Algorithm for calculating Born radii");
 948   EETYPE("gb-algorithm", ir->gb_algorithm, egb_names);
 949   CTYPE ("Frequency of calculating the Born radii inside rlist");
 950   ITYPE ("nstgbradii", ir->nstgbradii, 1);
 951   CTYPE ("Cutoff for Born radii calculation; the contribution from atoms");
 952   CTYPE ("between rlist and rgbradii is updated every nstlist steps");
 953   RTYPE ("rgbradii",  ir->rgbradii, 1.0);
 954   CTYPE ("Dielectric coefficient of the implicit solvent");
 955   RTYPE ("gb-epsilon-solvent",ir->gb_epsilon_solvent, 80.0);
 956   CTYPE ("Salt concentration in M for Generalized Born models");
 957   RTYPE ("gb-saltconc",  ir->gb_saltconc, 0.0);
 958   CTYPE ("Scaling factors used in the OBC GB model. Default values are OBC(II)");
 959   RTYPE ("gb-obc-alpha", ir->gb_obc_alpha, 1.0);
 960   RTYPE ("gb-obc-beta", ir->gb_obc_beta, 0.8);
 961   RTYPE ("gb-obc-gamma", ir->gb_obc_gamma, 4.85);
 962   RTYPE ("gb-dielectric-offset", ir->gb_dielectric_offset, 0.009);
 963   EETYPE("sa-algorithm", ir->sa_algorithm, esa_names);
 964   CTYPE ("Surface tension (kJ/mol/nm^2) for the SA (nonpolar surface) part of GBSA");
 965   CTYPE ("The value -1 will set default value for Still/HCT/OBC GB-models.");
 966   RTYPE ("sa-surface-tension", ir->sa_surface_tension, -1);
 967
 968   /* Coupling stuff */
 969   CCTYPE ("OPTIONS FOR WEAK COUPLING ALGORITHMS");
 970   CTYPE ("Temperature coupling");
 971   EETYPE("tcoupl",      ir->etc,        etcoupl_names);
 972   ITYPE ("nsttcouple", ir->nsttcouple,  -1);
 973   ITYPE("nh-chain-length",     ir->opts.nhchainlength, NHCHAINLENGTH);
 974   CTYPE ("Groups to couple separately");
 975   STYPE ("tc-grps",     tcgrps,         NULL);
 976   CTYPE ("Time constant (ps) and reference temperature (K)");
 977   STYPE ("tau-t",       tau_t,          NULL);
 978   STYPE ("ref-t",       ref_t,          NULL);
 979   CTYPE ("pressure coupling");
 980   EETYPE("pcoupl",      ir->epc,        epcoupl_names);
 981   EETYPE("pcoupltype",  ir->epct,       epcoupltype_names);
 982   ITYPE ("nstpcouple", ir->nstpcouple,  -1);
 983   CTYPE ("Time constant (ps), compressibility (1/bar) and reference P (bar)");
 984   RTYPE ("tau-p",       ir->tau_p,      1.0);
 985   STYPE ("compressibility",     dumstr[0],      NULL);
 986   STYPE ("ref-p",       dumstr[1],      NULL);
 987   CTYPE ("Scaling of reference coordinates, No, All or COM");
 988   EETYPE ("refcoord-scaling",ir->refcoord_scaling,erefscaling_names);
 989
 990   CTYPE ("Random seed for Andersen thermostat");
 991   ITYPE ("andersen-seed", ir->andersen_seed, 815131);
 992
 993   /* QMMM */
 994   CCTYPE ("OPTIONS FOR QMMM calculations");
 995   EETYPE("QMMM", ir->bQMMM, yesno_names);
 996   CTYPE ("Groups treated Quantum Mechanically");
 997   STYPE ("QMMM-grps",  QMMM,          NULL);
 998   CTYPE ("QM method");
 999   STYPE("QMmethod",     QMmethod, NULL);
1000   CTYPE ("QMMM scheme");
1001   EETYPE("QMMMscheme",  ir->QMMMscheme,    eQMMMscheme_names);
1002   CTYPE ("QM basisset");
1003   STYPE("QMbasis",      QMbasis, NULL);
1004   CTYPE ("QM charge");
1005   STYPE ("QMcharge",    QMcharge,NULL);
1006   CTYPE ("QM multiplicity");
1007   STYPE ("QMmult",      QMmult,NULL);
1008   CTYPE ("Surface Hopping");
1009   STYPE ("SH",          bSH, NULL);
1010   CTYPE ("CAS space options");
1011   STYPE ("CASorbitals",      CASorbitals,   NULL);
1012   STYPE ("CASelectrons",     CASelectrons,  NULL);
1013   STYPE ("SAon", SAon, NULL);
1014   STYPE ("SAoff",SAoff,NULL);
1015   STYPE ("SAsteps",  SAsteps, NULL);
1016   CTYPE ("Scale factor for MM charges");
1017   RTYPE ("MMChargeScaleFactor", ir->scalefactor, 1.0);
1018   CTYPE ("Optimization of QM subsystem");
1019   STYPE ("bOPT",          bOPT, NULL);
1020   STYPE ("bTS",          bTS, NULL);
1021
1022   /* Simulated annealing */
1023   CCTYPE("SIMULATED ANNEALING");
1024   CTYPE ("Type of annealing for each temperature group (no/single/periodic)");
1025   STYPE ("annealing",   anneal,      NULL);
1026   CTYPE ("Number of time points to use for specifying annealing in each group");
1027   STYPE ("annealing-npoints", anneal_npoints, NULL);
1028   CTYPE ("List of times at the annealing points for each group");
1029   STYPE ("annealing-time",       anneal_time,       NULL);
1030   CTYPE ("Temp. at each annealing point, for each group.");
1031   STYPE ("annealing-temp",  anneal_temp,  NULL);
1032
1033   /* Startup run */
1034   CCTYPE ("GENERATE VELOCITIES FOR STARTUP RUN");
1035   EETYPE("gen-vel",     opts->bGenVel,  yesno_names);
1036   RTYPE ("gen-temp",    opts->tempi,    300.0);
1037   ITYPE ("gen-seed",    opts->seed,     173529);
1038
1039   /* Shake stuff */
1040   CCTYPE ("OPTIONS FOR BONDS");
1041   EETYPE("constraints", opts->nshake,   constraints);
1042   CTYPE ("Type of constraint algorithm");
1043   EETYPE("constraint-algorithm",  ir->eConstrAlg, econstr_names);
1044   CTYPE ("Do not constrain the start configuration");
1045   EETYPE("continuation", ir->bContinuation, yesno_names);
1046   CTYPE ("Use successive overrelaxation to reduce the number of shake iterations");
1047   EETYPE("Shake-SOR", ir->bShakeSOR, yesno_names);
1048   CTYPE ("Relative tolerance of shake");
1049   RTYPE ("shake-tol", ir->shake_tol, 0.0001);
1050   CTYPE ("Highest order in the expansion of the constraint coupling matrix");
1051   ITYPE ("lincs-order", ir->nProjOrder, 4);
1052   CTYPE ("Number of iterations in the final step of LINCS. 1 is fine for");
1053   CTYPE ("normal simulations, but use 2 to conserve energy in NVE runs.");
1054   CTYPE ("For energy minimization with constraints it should be 4 to 8.");
1055   ITYPE ("lincs-iter", ir->nLincsIter, 1);
1056   CTYPE ("Lincs will write a warning to the stderr if in one step a bond");
1057   CTYPE ("rotates over more degrees than");
1058   RTYPE ("lincs-warnangle", ir->LincsWarnAngle, 30.0);
1059   CTYPE ("Convert harmonic bonds to morse potentials");
1060   EETYPE("morse",       opts->bMorse,yesno_names);
1061
1062   /* Energy group exclusions */
1063   CCTYPE ("ENERGY GROUP EXCLUSIONS");
1064   CTYPE ("Pairs of energy groups for which all non-bonded interactions are excluded");
1065   STYPE ("energygrp-excl", egpexcl,     NULL);
1066
1067   /* Walls */
1068   CCTYPE ("WALLS");
1069   CTYPE ("Number of walls, type, atom types, densities and box-z scale factor for Ewald");
1070   ITYPE ("nwall", ir->nwall, 0);
1071   EETYPE("wall-type",     ir->wall_type,   ewt_names);
1072   RTYPE ("wall-r-linpot", ir->wall_r_linpot, -1);
1073   STYPE ("wall-atomtype", wall_atomtype, NULL);
1074   STYPE ("wall-density",  wall_density,  NULL);
1075   RTYPE ("wall-ewald-zfac", ir->wall_ewald_zfac, 3);
1076
1077   /* COM pulling */
1078   CCTYPE("COM PULLING");
1079   CTYPE("Pull type: no, umbrella, constraint or constant-force");
1080   EETYPE("pull",          ir->ePull, epull_names);
1081   if (ir->ePull != epullNO) {
1082     snew(ir->pull,1);
1083     pull_grp = read_pullparams(&ninp,&inp,ir->pull,&opts->pull_start,wi);
1084   }
1085
1086   /* Enforced rotation */
1087   CCTYPE("ENFORCED ROTATION");
1088   CTYPE("Enforced rotation: No or Yes");
1089   EETYPE("rotation",       ir->bRot, yesno_names);
1090   if (ir->bRot) {
1091     snew(ir->rot,1);
1092     rot_grp = read_rotparams(&ninp,&inp,ir->rot,wi);
1093   }
1094
1095   /* Refinement */
1096   CCTYPE("NMR refinement stuff");
1097   CTYPE ("Distance restraints type: No, Simple or Ensemble");
1098   EETYPE("disre",       ir->eDisre,     edisre_names);
1099   CTYPE ("Force weighting of pairs in one distance restraint: Conservative or Equal");
1100   EETYPE("disre-weighting", ir->eDisreWeighting, edisreweighting_names);
1101   CTYPE ("Use sqrt of the time averaged times the instantaneous violation");
1102   EETYPE("disre-mixed", ir->bDisreMixed, yesno_names);
1103   RTYPE ("disre-fc",    ir->dr_fc,      1000.0);
1104   RTYPE ("disre-tau",   ir->dr_tau,     0.0);
1105   CTYPE ("Output frequency for pair distances to energy file");
1106   ITYPE ("nstdisreout", ir->nstdisreout, 100);
1107   CTYPE ("Orientation restraints: No or Yes");
1108   EETYPE("orire",       opts->bOrire,   yesno_names);
1109   CTYPE ("Orientation restraints force constant and tau for time averaging");
1110   RTYPE ("orire-fc",    ir->orires_fc,  0.0);
1111   RTYPE ("orire-tau",   ir->orires_tau, 0.0);
1112   STYPE ("orire-fitgrp",orirefitgrp,    NULL);
1113   CTYPE ("Output frequency for trace(SD) and S to energy file");
1114   ITYPE ("nstorireout", ir->nstorireout, 100);
1115   CTYPE ("Dihedral angle restraints: No or Yes");
1116   EETYPE("dihre",       opts->bDihre,   yesno_names);
1117   RTYPE ("dihre-fc",    ir->dihre_fc,   1000.0);
1118
1119   /* Free energy stuff */
1120   CCTYPE ("Free energy control stuff");
1121   EETYPE("free-energy", ir->efep, efep_names);
1122   RTYPE ("init-lambda", ir->init_lambda,0.0);
1123   RTYPE ("delta-lambda",ir->delta_lambda,0.0);
1124   STYPE ("foreign-lambda", foreign_lambda, NULL);
1125   RTYPE ("sc-alpha",ir->sc_alpha,0.0);
1126   ITYPE ("sc-power",ir->sc_power,0);
1127   RTYPE ("sc-sigma",ir->sc_sigma,0.3);
1128   ITYPE ("nstdhdl",     ir->nstdhdl, 10);
1129   EETYPE("separate-dhdl-file", ir->separate_dhdl_file,
1130                                separate_dhdl_file_names);
1131   EETYPE("dhdl-derivatives", ir->dhdl_derivatives, dhdl_derivatives_names);
1132   ITYPE ("dh-hist-size", ir->dh_hist_size, 0);
1133   RTYPE ("dh-hist-spacing", ir->dh_hist_spacing, 0.1);
1134   STYPE ("couple-moltype",  couple_moltype,  NULL);
1135   EETYPE("couple-lambda0", opts->couple_lam0, couple_lam);
1136   EETYPE("couple-lambda1", opts->couple_lam1, couple_lam);
1137   EETYPE("couple-intramol", opts->bCoupleIntra, yesno_names);
1138
1139   /* Non-equilibrium MD stuff */
1140   CCTYPE("Non-equilibrium MD stuff");
1141   STYPE ("acc-grps",    accgrps,        NULL);
1142   STYPE ("accelerate",  acc,            NULL);
1143   STYPE ("freezegrps",  freeze,         NULL);
1144   STYPE ("freezedim",   frdim,          NULL);
1145   RTYPE ("cos-acceleration", ir->cos_accel, 0);
1146   STYPE ("deform",      deform,         NULL);
1147
1148   /* Electric fields */
1149   CCTYPE("Electric fields");
1150   CTYPE ("Format is number of terms (int) and for all terms an amplitude (real)");
1151   CTYPE ("and a phase angle (real)");
1152   STYPE ("E-x",         efield_x,       NULL);
1153   STYPE ("E-xt",        efield_xt,      NULL);
1154   STYPE ("E-y",         efield_y,       NULL);
1155   STYPE ("E-yt",        efield_yt,      NULL);
1156   STYPE ("E-z",         efield_z,       NULL);
1157   STYPE ("E-zt",        efield_zt,      NULL);
1158
1159   /* User defined thingies */
1160   CCTYPE ("User defined thingies");
1161   STYPE ("user1-grps",  user1,          NULL);
1162   STYPE ("user2-grps",  user2,          NULL);
1163   ITYPE ("userint1",    ir->userint1,   0);
1164   ITYPE ("userint2",    ir->userint2,   0);
1165   ITYPE ("userint3",    ir->userint3,   0);
1166   ITYPE ("userint4",    ir->userint4,   0);
1167   RTYPE ("userreal1",   ir->userreal1,  0);
1168   RTYPE ("userreal2",   ir->userreal2,  0);
1169   RTYPE ("userreal3",   ir->userreal3,  0);
1170   RTYPE ("userreal4",   ir->userreal4,  0);
1171 #undef CTYPE
1172
1173   write_inpfile(mdparout,ninp,inp,FALSE,wi);
1174   for (i=0; (i<ninp); i++) {
1175     sfree(inp[i].name);
1176     sfree(inp[i].value);
1177   }
1178   sfree(inp);
1179
1180   /* Process options if necessary */
1181   for(m=0; m<2; m++) {
1182     for(i=0; i<2*DIM; i++)
1183       dumdub[m][i]=0.0;
1184     if(ir->epc) {
1185       switch (ir->epct) {
1186       case epctISOTROPIC:
1187         if (sscanf(dumstr[m],"%lf",&(dumdub[m][XX]))!=1) {
1188         warning_error(wi,"Pressure coupling not enough values (I need 1)");
1189         }
1190         dumdub[m][YY]=dumdub[m][ZZ]=dumdub[m][XX];
1191         break;
1192       case epctSEMIISOTROPIC:
1193       case epctSURFACETENSION:
1194         if (sscanf(dumstr[m],"%lf%lf",
1195                    &(dumdub[m][XX]),&(dumdub[m][ZZ]))!=2) {
1196         warning_error(wi,"Pressure coupling not enough values (I need 2)");
1197         }
1198         dumdub[m][YY]=dumdub[m][XX];
1199         break;
1200       case epctANISOTROPIC:
1201         if (sscanf(dumstr[m],"%lf%lf%lf%lf%lf%lf",
1202                    &(dumdub[m][XX]),&(dumdub[m][YY]),&(dumdub[m][ZZ]),
1203                    &(dumdub[m][3]),&(dumdub[m][4]),&(dumdub[m][5]))!=6) {
1204         warning_error(wi,"Pressure coupling not enough values (I need 6)");
1205         }
1206         break;
1207       default:
1208         gmx_fatal(FARGS,"Pressure coupling type %s not implemented yet",
1209                     epcoupltype_names[ir->epct]);
1210       }
1211     }
1212   }
1213   clear_mat(ir->ref_p);
1214   clear_mat(ir->compress);
1215   for(i=0; i<DIM; i++) {
1216     ir->ref_p[i][i]    = dumdub[1][i];
1217     ir->compress[i][i] = dumdub[0][i];
1218   }
1219   if (ir->epct == epctANISOTROPIC) {
1220     ir->ref_p[XX][YY] = dumdub[1][3];
1221     ir->ref_p[XX][ZZ] = dumdub[1][4];
1222     ir->ref_p[YY][ZZ] = dumdub[1][5];
1223     if (ir->ref_p[XX][YY]!=0 && ir->ref_p[XX][ZZ]!=0 && ir->ref_p[YY][ZZ]!=0) {
1224       warning(wi,"All off-diagonal reference pressures are non-zero. Are you sure you want to apply a threefold shear stress?\n");
1225     }
1226     ir->compress[XX][YY] = dumdub[0][3];
1227     ir->compress[XX][ZZ] = dumdub[0][4];
1228     ir->compress[YY][ZZ] = dumdub[0][5];
1229     for(i=0; i<DIM; i++) {
1230       for(m=0; m<i; m++) {
1231         ir->ref_p[i][m] = ir->ref_p[m][i];
1232         ir->compress[i][m] = ir->compress[m][i];
1233       }
1234     }
1235   }
1236
1237   if (ir->comm_mode == ecmNO)
1238     ir->nstcomm = 0;
1239
1240   opts->couple_moltype = NULL;
1241   if (strlen(couple_moltype) > 0) {
1242     if (ir->efep != efepNO) {
1243       opts->couple_moltype = strdup(couple_moltype);
1244       if (opts->couple_lam0 == opts->couple_lam1)
1245         warning(wi,"The lambda=0 and lambda=1 states for coupling are identical");
1246       if (ir->eI == eiMD && (opts->couple_lam0 == ecouplamNONE ||
1247                              opts->couple_lam1 == ecouplamNONE)) {
1248         warning(wi,"For proper sampling of the (nearly) decoupled state, stochastic dynamics should be used");
1249       }
1250     } else {
1251       warning(wi,"Can not couple a molecule with free-energy = no");
1252     }
1253   }
1254
1255   do_wall_params(ir,wall_atomtype,wall_density,opts);
1256
1257   if (opts->bOrire && str_nelem(orirefitgrp,MAXPTR,NULL)!=1) {
1258       warning_error(wi,"ERROR: Need one orientation restraint fit group\n");
1259   }
1260
1261   clear_mat(ir->deform);
1262   for(i=0; i<6; i++)
1263     dumdub[0][i] = 0;
1264   m = sscanf(deform,"%lf %lf %lf %lf %lf %lf",
1265              &(dumdub[0][0]),&(dumdub[0][1]),&(dumdub[0][2]),
1266              &(dumdub[0][3]),&(dumdub[0][4]),&(dumdub[0][5]));
1267   for(i=0; i<3; i++)
1268     ir->deform[i][i] = dumdub[0][i];
1269   ir->deform[YY][XX] = dumdub[0][3];
1270   ir->deform[ZZ][XX] = dumdub[0][4];
1271   ir->deform[ZZ][YY] = dumdub[0][5];
1272   if (ir->epc != epcNO) {
1273     for(i=0; i<3; i++)
1274       for(j=0; j<=i; j++)
1275         if (ir->deform[i][j]!=0 && ir->compress[i][j]!=0) {
1276         warning_error(wi,"A box element has deform set and compressibility > 0");
1277         }
1278     for(i=0; i<3; i++)
1279       for(j=0; j<i; j++)
1280         if (ir->deform[i][j]!=0) {
1281           for(m=j; m<DIM; m++)
1282             if (ir->compress[m][j]!=0) {
1283               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.");
1284               warning(wi,warn_buf);
1285             }
1286         }
1287   }
1288
1289   if (ir->efep != efepNO) {
1290     parse_n_double(foreign_lambda,&ir->n_flambda,&ir->flambda);
1291     if (ir->n_flambda > 0 && ir->rlist < max(ir->rvdw,ir->rcoulomb)) {
1292       warning_note(wi,"For foreign lambda free energy differences it is assumed that the soft-core interactions have no effect beyond the neighborlist cut-off");
1293     }
1294   } else {
1295     ir->n_flambda = 0;
1296   }
1297
1298   sfree(dumstr[0]);
1299   sfree(dumstr[1]);
1300 }
1301
1302 static int search_QMstring(char *s,int ng,const char *gn[])
1303 {
1304   /* same as normal search_string, but this one searches QM strings */
1305   int i;
1306
1307   for(i=0; (i<ng); i++)
1308     if (gmx_strcasecmp(s,gn[i]) == 0)
1309       return i;
1310
1311   gmx_fatal(FARGS,"this QM method or basisset (%s) is not implemented\n!",s);
1312
1313   return -1;
1314
1315 } /* search_QMstring */
1316
1317
1318 int search_string(char *s,int ng,char *gn[])
1319 {
1320   int i;
1321
1322   for(i=0; (i<ng); i++)
1323   {
1324     if (gmx_strcasecmp(s,gn[i]) == 0)
1325     {
1326       return i;
1327     }
1328   }
1329
1330   gmx_fatal(FARGS,"Group %s not found in index file.\nGroup names must match either [moleculetype] names\nor custom index group names,in which case you\nmust supply an index file to the '-n' option of grompp.",s);
1331
1332   return -1;
1333 }
1334
1335 static gmx_bool do_numbering(int natoms,gmx_groups_t *groups,int ng,char *ptrs[],
1336                          t_blocka *block,char *gnames[],
1337                          int gtype,int restnm,
1338                          int grptp,gmx_bool bVerbose,
1339                          warninp_t wi)
1340 {
1341     unsigned short *cbuf;
1342     t_grps *grps=&(groups->grps[gtype]);
1343     int    i,j,gid,aj,ognr,ntot=0;
1344     const char *title;
1345     gmx_bool   bRest;
1346     char   warn_buf[STRLEN];
1347
1348     if (debug)
1349     {
1350         fprintf(debug,"Starting numbering %d groups of type %d\n",ng,gtype);
1351     }
1352
1353     title = gtypes[gtype];
1354
1355     snew(cbuf,natoms);
1356     /* Mark all id's as not set */
1357     for(i=0; (i<natoms); i++)
1358     {
1359         cbuf[i] = NOGID;
1360     }
1361
1362     snew(grps->nm_ind,ng+1); /* +1 for possible rest group */
1363     for(i=0; (i<ng); i++)
1364     {
1365         /* Lookup the group name in the block structure */
1366         gid = search_string(ptrs[i],block->nr,gnames);
1367         if ((grptp != egrptpONE) || (i == 0))
1368         {
1369             grps->nm_ind[grps->nr++]=gid;
1370         }
1371         if (debug)
1372         {
1373             fprintf(debug,"Found gid %d for group %s\n",gid,ptrs[i]);
1374         }
1375
1376         /* Now go over the atoms in the group */
1377         for(j=block->index[gid]; (j<block->index[gid+1]); j++)
1378         {
1379
1380             aj=block->a[j];
1381
1382             /* Range checking */
1383             if ((aj < 0) || (aj >= natoms))
1384             {
1385                 gmx_fatal(FARGS,"Invalid atom number %d in indexfile",aj);
1386             }
1387             /* Lookup up the old group number */
1388             ognr = cbuf[aj];
1389             if (ognr != NOGID)
1390             {
1391                 gmx_fatal(FARGS,"Atom %d in multiple %s groups (%d and %d)",
1392                           aj+1,title,ognr+1,i+1);
1393             }
1394             else
1395             {
1396                 /* Store the group number in buffer */
1397                 if (grptp == egrptpONE)
1398                 {
1399                     cbuf[aj] = 0;
1400                 }
1401                 else
1402                 {
1403                     cbuf[aj] = i;
1404                 }
1405                 ntot++;
1406             }
1407         }
1408     }
1409
1410     /* Now check whether we have done all atoms */
1411     bRest = FALSE;
1412     if (ntot != natoms)
1413     {
1414         if (grptp == egrptpALL)
1415         {
1416             gmx_fatal(FARGS,"%d atoms are not part of any of the %s groups",
1417                       natoms-ntot,title);
1418         }
1419         else if (grptp == egrptpPART)
1420         {
1421             sprintf(warn_buf,"%d atoms are not part of any of the %s groups",
1422                     natoms-ntot,title);
1423             warning_note(wi,warn_buf);
1424         }
1425         /* Assign all atoms currently unassigned to a rest group */
1426         for(j=0; (j<natoms); j++)
1427         {
1428             if (cbuf[j] == NOGID)
1429             {
1430                 cbuf[j] = grps->nr;
1431                 bRest = TRUE;
1432             }
1433         }
1434         if (grptp != egrptpPART)
1435         {
1436             if (bVerbose)
1437             {
1438                 fprintf(stderr,
1439                         "Making dummy/rest group for %s containing %d elements\n",
1440                         title,natoms-ntot);
1441             }
1442             /* Add group name "rest" */
1443             grps->nm_ind[grps->nr] = restnm;
1444
1445             /* Assign the rest name to all atoms not currently assigned to a group */
1446             for(j=0; (j<natoms); j++)
1447             {
1448                 if (cbuf[j] == NOGID)
1449                 {
1450                     cbuf[j] = grps->nr;
1451                 }
1452             }
1453             grps->nr++;
1454         }
1455     }
1456
1457     if (grps->nr == 1)
1458     {
1459         groups->ngrpnr[gtype] = 0;
1460         groups->grpnr[gtype]  = NULL;
1461     }
1462     else
1463     {
1464         groups->ngrpnr[gtype] = natoms;
1465         snew(groups->grpnr[gtype],natoms);
1466         for(j=0; (j<natoms); j++)
1467         {
1468             groups->grpnr[gtype][j] = cbuf[j];
1469         }
1470     }
1471
1472     sfree(cbuf);
1473
1474     return (bRest && grptp == egrptpPART);
1475 }
1476
1477 static void calc_nrdf(gmx_mtop_t *mtop,t_inputrec *ir,char **gnames)
1478 {
1479   t_grpopts *opts;
1480   gmx_groups_t *groups;
1481   t_pull  *pull;
1482   int     natoms,ai,aj,i,j,d,g,imin,jmin,nc;
1483   t_iatom *ia;
1484   int     *nrdf2,*na_vcm,na_tot;
1485   double  *nrdf_tc,*nrdf_vcm,nrdf_uc,n_sub=0;
1486   gmx_mtop_atomloop_all_t aloop;
1487   t_atom  *atom;
1488   int     mb,mol,ftype,as;
1489   gmx_molblock_t *molb;
1490   gmx_moltype_t *molt;
1491
1492   /* Calculate nrdf.
1493    * First calc 3xnr-atoms for each group
1494    * then subtract half a degree of freedom for each constraint
1495    *
1496    * Only atoms and nuclei contribute to the degrees of freedom...
1497    */
1498
1499   opts = &ir->opts;
1500
1501   groups = &mtop->groups;
1502   natoms = mtop->natoms;
1503
1504   /* Allocate one more for a possible rest group */
1505   /* We need to sum degrees of freedom into doubles,
1506    * since floats give too low nrdf's above 3 million atoms.
1507    */
1508   snew(nrdf_tc,groups->grps[egcTC].nr+1);
1509   snew(nrdf_vcm,groups->grps[egcVCM].nr+1);
1510   snew(na_vcm,groups->grps[egcVCM].nr+1);
1511
1512   for(i=0; i<groups->grps[egcTC].nr; i++)
1513     nrdf_tc[i] = 0;
1514   for(i=0; i<groups->grps[egcVCM].nr+1; i++)
1515     nrdf_vcm[i] = 0;
1516
1517   snew(nrdf2,natoms);
1518   aloop = gmx_mtop_atomloop_all_init(mtop);
1519   while (gmx_mtop_atomloop_all_next(aloop,&i,&atom)) {
1520     nrdf2[i] = 0;
1521     if (atom->ptype == eptAtom || atom->ptype == eptNucleus) {
1522       g = ggrpnr(groups,egcFREEZE,i);
1523       /* Double count nrdf for particle i */
1524       for(d=0; d<DIM; d++) {
1525         if (opts->nFreeze[g][d] == 0) {
1526           nrdf2[i] += 2;
1527         }
1528       }
1529       nrdf_tc [ggrpnr(groups,egcTC ,i)] += 0.5*nrdf2[i];
1530       nrdf_vcm[ggrpnr(groups,egcVCM,i)] += 0.5*nrdf2[i];
1531     }
1532   }
1533
1534   as = 0;
1535   for(mb=0; mb<mtop->nmolblock; mb++) {
1536     molb = &mtop->molblock[mb];
1537     molt = &mtop->moltype[molb->type];
1538     atom = molt->atoms.atom;
1539     for(mol=0; mol<molb->nmol; mol++) {
1540       for (ftype=F_CONSTR; ftype<=F_CONSTRNC; ftype++) {
1541         ia = molt->ilist[ftype].iatoms;
1542         for(i=0; i<molt->ilist[ftype].nr; ) {
1543           /* Subtract degrees of freedom for the constraints,
1544            * if the particles still have degrees of freedom left.
1545            * If one of the particles is a vsite or a shell, then all
1546            * constraint motion will go there, but since they do not
1547            * contribute to the constraints the degrees of freedom do not
1548            * change.
1549            */
1550           ai = as + ia[1];
1551           aj = as + ia[2];
1552           if (((atom[ia[1]].ptype == eptNucleus) ||
1553                (atom[ia[1]].ptype == eptAtom)) &&
1554               ((atom[ia[2]].ptype == eptNucleus) ||
1555                (atom[ia[2]].ptype == eptAtom))) {
1556             if (nrdf2[ai] > 0)
1557               jmin = 1;
1558             else
1559               jmin = 2;
1560             if (nrdf2[aj] > 0)
1561               imin = 1;
1562             else
1563               imin = 2;
1564             imin = min(imin,nrdf2[ai]);
1565             jmin = min(jmin,nrdf2[aj]);
1566             nrdf2[ai] -= imin;
1567             nrdf2[aj] -= jmin;
1568             nrdf_tc [ggrpnr(groups,egcTC ,ai)] -= 0.5*imin;
1569             nrdf_tc [ggrpnr(groups,egcTC ,aj)] -= 0.5*jmin;
1570             nrdf_vcm[ggrpnr(groups,egcVCM,ai)] -= 0.5*imin;
1571             nrdf_vcm[ggrpnr(groups,egcVCM,aj)] -= 0.5*jmin;
1572           }
1573           ia += interaction_function[ftype].nratoms+1;
1574           i  += interaction_function[ftype].nratoms+1;
1575         }
1576       }
1577       ia = molt->ilist[F_SETTLE].iatoms;
1578       for(i=0; i<molt->ilist[F_SETTLE].nr; ) {
1579         /* Subtract 1 dof from every atom in the SETTLE */
1580         for(ai=as+ia[1]; ai<as+ia[1]+3; ai++) {
1581           imin = min(2,nrdf2[ai]);
1582           nrdf2[ai] -= imin;
1583           nrdf_tc [ggrpnr(groups,egcTC ,ai)] -= 0.5*imin;
1584           nrdf_vcm[ggrpnr(groups,egcVCM,ai)] -= 0.5*imin;
1585         }
1586         ia += 2;
1587         i  += 2;
1588       }
1589       as += molt->atoms.nr;
1590     }
1591   }
1592
1593   if (ir->ePull == epullCONSTRAINT) {
1594     /* Correct nrdf for the COM constraints.
1595      * We correct using the TC and VCM group of the first atom
1596      * in the reference and pull group. If atoms in one pull group
1597      * belong to different TC or VCM groups it is anyhow difficult
1598      * to determine the optimal nrdf assignment.
1599      */
1600     pull = ir->pull;
1601     if (pull->eGeom == epullgPOS) {
1602       nc = 0;
1603       for(i=0; i<DIM; i++) {
1604         if (pull->dim[i])
1605           nc++;
1606       }
1607     } else {
1608       nc = 1;
1609     }
1610     for(i=0; i<pull->ngrp; i++) {
1611       imin = 2*nc;
1612       if (pull->grp[0].nat > 0) {
1613         /* Subtract 1/2 dof from the reference group */
1614         ai = pull->grp[0].ind[0];
1615         if (nrdf_tc[ggrpnr(groups,egcTC,ai)] > 1) {
1616           nrdf_tc [ggrpnr(groups,egcTC ,ai)] -= 0.5;
1617           nrdf_vcm[ggrpnr(groups,egcVCM,ai)] -= 0.5;
1618           imin--;
1619         }
1620       }
1621       /* Subtract 1/2 dof from the pulled group */
1622       ai = pull->grp[1+i].ind[0];
1623       nrdf_tc [ggrpnr(groups,egcTC ,ai)] -= 0.5*imin;
1624       nrdf_vcm[ggrpnr(groups,egcVCM,ai)] -= 0.5*imin;
1625       if (nrdf_tc[ggrpnr(groups,egcTC,ai)] < 0)
1626         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)]]);
1627     }
1628   }
1629
1630   if (ir->nstcomm != 0) {
1631     /* Subtract 3 from the number of degrees of freedom in each vcm group
1632      * when com translation is removed and 6 when rotation is removed
1633      * as well.
1634      */
1635     switch (ir->comm_mode) {
1636     case ecmLINEAR:
1637       n_sub = ndof_com(ir);
1638       break;
1639     case ecmANGULAR:
1640       n_sub = 6;
1641       break;
1642     default:
1643       n_sub = 0;
1644       gmx_incons("Checking comm_mode");
1645     }
1646
1647     for(i=0; i<groups->grps[egcTC].nr; i++) {
1648       /* Count the number of atoms of TC group i for every VCM group */
1649       for(j=0; j<groups->grps[egcVCM].nr+1; j++)
1650         na_vcm[j] = 0;
1651       na_tot = 0;
1652       for(ai=0; ai<natoms; ai++)
1653         if (ggrpnr(groups,egcTC,ai) == i) {
1654           na_vcm[ggrpnr(groups,egcVCM,ai)]++;
1655           na_tot++;
1656         }
1657       /* Correct for VCM removal according to the fraction of each VCM
1658        * group present in this TC group.
1659        */
1660       nrdf_uc = nrdf_tc[i];
1661       if (debug) {
1662         fprintf(debug,"T-group[%d] nrdf_uc = %g, n_sub = %g\n",
1663                 i,nrdf_uc,n_sub);
1664       }
1665       nrdf_tc[i] = 0;
1666       for(j=0; j<groups->grps[egcVCM].nr+1; j++) {
1667         if (nrdf_vcm[j] > n_sub) {
1668           nrdf_tc[i] += nrdf_uc*((double)na_vcm[j]/(double)na_tot)*
1669             (nrdf_vcm[j] - n_sub)/nrdf_vcm[j];
1670         }
1671         if (debug) {
1672           fprintf(debug,"  nrdf_vcm[%d] = %g, nrdf = %g\n",
1673                   j,nrdf_vcm[j],nrdf_tc[i]);
1674         }
1675       }
1676     }
1677   }
1678   for(i=0; (i<groups->grps[egcTC].nr); i++) {
1679     opts->nrdf[i] = nrdf_tc[i];
1680     if (opts->nrdf[i] < 0)
1681       opts->nrdf[i] = 0;
1682     fprintf(stderr,
1683             "Number of degrees of freedom in T-Coupling group %s is %.2f\n",
1684             gnames[groups->grps[egcTC].nm_ind[i]],opts->nrdf[i]);
1685   }
1686
1687   sfree(nrdf2);
1688   sfree(nrdf_tc);
1689   sfree(nrdf_vcm);
1690   sfree(na_vcm);
1691 }
1692
1693 static void decode_cos(char *s,t_cosines *cosine,gmx_bool bTime)
1694 {
1695   char   *t;
1696   char   format[STRLEN],f1[STRLEN];
1697   double a,phi;
1698   int    i;
1699
1700   t=strdup(s);
1701   trim(t);
1702
1703   cosine->n=0;
1704   cosine->a=NULL;
1705   cosine->phi=NULL;
1706   if (strlen(t)) {
1707     sscanf(t,"%d",&(cosine->n));
1708     if (cosine->n <= 0) {
1709       cosine->n=0;
1710     } else {
1711       snew(cosine->a,cosine->n);
1712       snew(cosine->phi,cosine->n);
1713
1714       sprintf(format,"%%*d");
1715       for(i=0; (i<cosine->n); i++) {
1716         strcpy(f1,format);
1717         strcat(f1,"%lf%lf");
1718         if (sscanf(t,f1,&a,&phi) < 2)
1719           gmx_fatal(FARGS,"Invalid input for electric field shift: '%s'",t);
1720         cosine->a[i]=a;
1721         cosine->phi[i]=phi;
1722         strcat(format,"%*lf%*lf");
1723       }
1724     }
1725   }
1726   sfree(t);
1727 }
1728
1729 static gmx_bool do_egp_flag(t_inputrec *ir,gmx_groups_t *groups,
1730                         const char *option,const char *val,int flag)
1731 {
1732   /* The maximum number of energy group pairs would be MAXPTR*(MAXPTR+1)/2.
1733    * But since this is much larger than STRLEN, such a line can not be parsed.
1734    * The real maximum is the number of names that fit in a string: STRLEN/2.
1735    */
1736 #define EGP_MAX (STRLEN/2)
1737   int  nelem,i,j,k,nr;
1738   char *names[EGP_MAX];
1739   char ***gnames;
1740   gmx_bool bSet;
1741
1742   gnames = groups->grpname;
1743
1744   nelem = str_nelem(val,EGP_MAX,names);
1745   if (nelem % 2 != 0)
1746     gmx_fatal(FARGS,"The number of groups for %s is odd",option);
1747   nr = groups->grps[egcENER].nr;
1748   bSet = FALSE;
1749   for(i=0; i<nelem/2; i++) {
1750     j = 0;
1751     while ((j < nr) &&
1752            gmx_strcasecmp(names[2*i],*(gnames[groups->grps[egcENER].nm_ind[j]])))
1753       j++;
1754     if (j == nr)
1755       gmx_fatal(FARGS,"%s in %s is not an energy group\n",
1756                   names[2*i],option);
1757     k = 0;
1758     while ((k < nr) &&
1759            gmx_strcasecmp(names[2*i+1],*(gnames[groups->grps[egcENER].nm_ind[k]])))
1760       k++;
1761     if (k==nr)
1762       gmx_fatal(FARGS,"%s in %s is not an energy group\n",
1763               names[2*i+1],option);
1764     if ((j < nr) && (k < nr)) {
1765       ir->opts.egp_flags[nr*j+k] |= flag;
1766       ir->opts.egp_flags[nr*k+j] |= flag;
1767       bSet = TRUE;
1768     }
1769   }
1770
1771   return bSet;
1772 }
1773
1774 void do_index(const char* mdparin, const char *ndx,
1775               gmx_mtop_t *mtop,
1776               gmx_bool bVerbose,
1777               t_inputrec *ir,rvec *v,
1778               warninp_t wi)
1779 {
1780   t_blocka *grps;
1781   gmx_groups_t *groups;
1782   int     natoms;
1783   t_symtab *symtab;
1784   t_atoms atoms_all;
1785   char    warnbuf[STRLEN],**gnames;
1786   int     nr,ntcg,ntau_t,nref_t,nacc,nofg,nSA,nSA_points,nSA_time,nSA_temp;
1787   real    tau_min;
1788   int     nstcmin;
1789   int     nacg,nfreeze,nfrdim,nenergy,nvcm,nuser;
1790   char    *ptr1[MAXPTR],*ptr2[MAXPTR],*ptr3[MAXPTR];
1791   int     i,j,k,restnm;
1792   real    SAtime;
1793   gmx_bool    bExcl,bTable,bSetTCpar,bAnneal,bRest;
1794   int     nQMmethod,nQMbasis,nQMcharge,nQMmult,nbSH,nCASorb,nCASelec,
1795     nSAon,nSAoff,nSAsteps,nQMg,nbOPT,nbTS;
1796   char    warn_buf[STRLEN];
1797
1798   if (bVerbose)
1799     fprintf(stderr,"processing index file...\n");
1800   debug_gmx();
1801   if (ndx == NULL) {
1802     snew(grps,1);
1803     snew(grps->index,1);
1804     snew(gnames,1);
1805     atoms_all = gmx_mtop_global_atoms(mtop);
1806     analyse(&atoms_all,grps,&gnames,FALSE,TRUE);
1807     free_t_atoms(&atoms_all,FALSE);
1808   } else {
1809     grps = init_index(ndx,&gnames);
1810   }
1811
1812   groups = &mtop->groups;
1813   natoms = mtop->natoms;
1814   symtab = &mtop->symtab;
1815
1816   snew(groups->grpname,grps->nr+1);
1817
1818   for(i=0; (i<grps->nr); i++) {
1819     groups->grpname[i] = put_symtab(symtab,gnames[i]);
1820   }
1821   groups->grpname[i] = put_symtab(symtab,"rest");
1822   restnm=i;
1823   srenew(gnames,grps->nr+1);
1824   gnames[restnm] = *(groups->grpname[i]);
1825   groups->ngrpname = grps->nr+1;
1826
1827   set_warning_line(wi,mdparin,-1);
1828
1829   ntau_t = str_nelem(tau_t,MAXPTR,ptr1);
1830   nref_t = str_nelem(ref_t,MAXPTR,ptr2);
1831   ntcg   = str_nelem(tcgrps,MAXPTR,ptr3);
1832   if ((ntau_t != ntcg) || (nref_t != ntcg)) {
1833     gmx_fatal(FARGS,"Invalid T coupling input: %d groups, %d ref-t values and "
1834                 "%d tau-t values",ntcg,nref_t,ntau_t);
1835   }
1836
1837   bSetTCpar = (ir->etc || EI_SD(ir->eI) || ir->eI==eiBD || EI_TPI(ir->eI));
1838   do_numbering(natoms,groups,ntcg,ptr3,grps,gnames,egcTC,
1839                restnm,bSetTCpar ? egrptpALL : egrptpALL_GENREST,bVerbose,wi);
1840   nr = groups->grps[egcTC].nr;
1841   ir->opts.ngtc = nr;
1842   snew(ir->opts.nrdf,nr);
1843   snew(ir->opts.tau_t,nr);
1844   snew(ir->opts.ref_t,nr);
1845   if (ir->eI==eiBD && ir->bd_fric==0) {
1846     fprintf(stderr,"bd-fric=0, so tau-t will be used as the inverse friction constant(s)\n");
1847   }
1848
1849   if (bSetTCpar)
1850   {
1851       if (nr != nref_t)
1852       {
1853           gmx_fatal(FARGS,"Not enough ref-t and tau-t values!");
1854       }
1855
1856       tau_min = 1e20;
1857       for(i=0; (i<nr); i++)
1858       {
1859           ir->opts.tau_t[i] = strtod(ptr1[i],NULL);
1860           if ((ir->eI == eiBD || ir->eI == eiSD2) && ir->opts.tau_t[i] <= 0)
1861           {
1862               sprintf(warn_buf,"With integrator %s tau-t should be larger than 0",ei_names[ir->eI]);
1863               warning_error(wi,warn_buf);
1864           }
1865           if ((ir->etc == etcVRESCALE && ir->opts.tau_t[i] >= 0) ||
1866               (ir->etc != etcVRESCALE && ir->opts.tau_t[i] >  0))
1867           {
1868               tau_min = min(tau_min,ir->opts.tau_t[i]);
1869           }
1870       }
1871       if (ir->etc != etcNO && ir->nsttcouple == -1)
1872       {
1873             ir->nsttcouple = ir_optimal_nsttcouple(ir);
1874       }
1875       if (EI_VV(ir->eI))
1876       {
1877           if ((ir->epc==epcMTTK) && (ir->etc>etcNO))
1878           {
1879               int mincouple;
1880               mincouple = ir->nsttcouple;
1881               if (ir->nstpcouple < mincouple)
1882               {
1883                   mincouple = ir->nstpcouple;
1884               }
1885               ir->nstpcouple = mincouple;
1886               ir->nsttcouple = mincouple;
1887               warning_note(wi,"for current Trotter decomposition methods with vv, nsttcouple and nstpcouple must be equal.  Both have been reset to min(nsttcouple,nstpcouple)");
1888           }
1889       }
1890       nstcmin = tcouple_min_integration_steps(ir->etc);
1891       if (nstcmin > 1)
1892       {
1893           if (tau_min/(ir->delta_t*ir->nsttcouple) < nstcmin)
1894           {
1895               sprintf(warn_buf,"For proper integration of the %s thermostat, tau-t (%g) should be at least %d times larger than nsttcouple*dt (%g)",
1896                       ETCOUPLTYPE(ir->etc),
1897                       tau_min,nstcmin,
1898                       ir->nsttcouple*ir->delta_t);
1899               warning(wi,warn_buf);
1900           }
1901       }
1902       for(i=0; (i<nr); i++)
1903       {
1904           ir->opts.ref_t[i] = strtod(ptr2[i],NULL);
1905           if (ir->opts.ref_t[i] < 0)
1906           {
1907               gmx_fatal(FARGS,"ref-t for group %d negative",i);
1908           }
1909       }
1910   }
1911
1912   /* Simulated annealing for each group. There are nr groups */
1913   nSA = str_nelem(anneal,MAXPTR,ptr1);
1914   if (nSA == 1 && (ptr1[0][0]=='n' || ptr1[0][0]=='N'))
1915      nSA = 0;
1916   if(nSA>0 && nSA != nr)
1917     gmx_fatal(FARGS,"Not enough annealing values: %d (for %d groups)\n",nSA,nr);
1918   else {
1919     snew(ir->opts.annealing,nr);
1920     snew(ir->opts.anneal_npoints,nr);
1921     snew(ir->opts.anneal_time,nr);
1922     snew(ir->opts.anneal_temp,nr);
1923     for(i=0;i<nr;i++) {
1924       ir->opts.annealing[i]=eannNO;
1925       ir->opts.anneal_npoints[i]=0;
1926       ir->opts.anneal_time[i]=NULL;
1927       ir->opts.anneal_temp[i]=NULL;
1928     }
1929     if (nSA > 0) {
1930       bAnneal=FALSE;
1931       for(i=0;i<nr;i++) {
1932         if(ptr1[i][0]=='n' || ptr1[i][0]=='N') {
1933           ir->opts.annealing[i]=eannNO;
1934         } else if(ptr1[i][0]=='s'|| ptr1[i][0]=='S') {
1935           ir->opts.annealing[i]=eannSINGLE;
1936           bAnneal=TRUE;
1937         } else if(ptr1[i][0]=='p'|| ptr1[i][0]=='P') {
1938           ir->opts.annealing[i]=eannPERIODIC;
1939           bAnneal=TRUE;
1940         }
1941       }
1942       if(bAnneal) {
1943         /* Read the other fields too */
1944         nSA_points = str_nelem(anneal_npoints,MAXPTR,ptr1);
1945         if(nSA_points!=nSA)
1946           gmx_fatal(FARGS,"Found %d annealing-npoints values for %d groups\n",nSA_points,nSA);
1947         for(k=0,i=0;i<nr;i++) {
1948           ir->opts.anneal_npoints[i]=strtol(ptr1[i],NULL,10);
1949           if(ir->opts.anneal_npoints[i]==1)
1950             gmx_fatal(FARGS,"Please specify at least a start and an end point for annealing\n");
1951           snew(ir->opts.anneal_time[i],ir->opts.anneal_npoints[i]);
1952           snew(ir->opts.anneal_temp[i],ir->opts.anneal_npoints[i]);
1953           k += ir->opts.anneal_npoints[i];
1954         }
1955
1956         nSA_time = str_nelem(anneal_time,MAXPTR,ptr1);
1957         if(nSA_time!=k)
1958           gmx_fatal(FARGS,"Found %d annealing-time values, wanter %d\n",nSA_time,k);
1959         nSA_temp = str_nelem(anneal_temp,MAXPTR,ptr2);
1960         if(nSA_temp!=k)
1961           gmx_fatal(FARGS,"Found %d annealing-temp values, wanted %d\n",nSA_temp,k);
1962
1963         for(i=0,k=0;i<nr;i++) {
1964
1965           for(j=0;j<ir->opts.anneal_npoints[i];j++) {
1966             ir->opts.anneal_time[i][j]=strtod(ptr1[k],NULL);
1967             ir->opts.anneal_temp[i][j]=strtod(ptr2[k],NULL);
1968             if(j==0) {
1969               if(ir->opts.anneal_time[i][0] > (ir->init_t+GMX_REAL_EPS))
1970                 gmx_fatal(FARGS,"First time point for annealing > init_t.\n");
1971             } else {
1972               /* j>0 */
1973               if(ir->opts.anneal_time[i][j]<ir->opts.anneal_time[i][j-1])
1974                 gmx_fatal(FARGS,"Annealing timepoints out of order: t=%f comes after t=%f\n",
1975                             ir->opts.anneal_time[i][j],ir->opts.anneal_time[i][j-1]);
1976             }
1977             if(ir->opts.anneal_temp[i][j]<0)
1978               gmx_fatal(FARGS,"Found negative temperature in annealing: %f\n",ir->opts.anneal_temp[i][j]);
1979             k++;
1980           }
1981         }
1982         /* Print out some summary information, to make sure we got it right */
1983         for(i=0,k=0;i<nr;i++) {
1984           if(ir->opts.annealing[i]!=eannNO) {
1985             j = groups->grps[egcTC].nm_ind[i];
1986             fprintf(stderr,"Simulated annealing for group %s: %s, %d timepoints\n",
1987                     *(groups->grpname[j]),eann_names[ir->opts.annealing[i]],
1988                     ir->opts.anneal_npoints[i]);
1989             fprintf(stderr,"Time (ps)   Temperature (K)\n");
1990             /* All terms except the last one */
1991             for(j=0;j<(ir->opts.anneal_npoints[i]-1);j++)
1992                 fprintf(stderr,"%9.1f      %5.1f\n",ir->opts.anneal_time[i][j],ir->opts.anneal_temp[i][j]);
1993
1994             /* Finally the last one */
1995             j = ir->opts.anneal_npoints[i]-1;
1996             if(ir->opts.annealing[i]==eannSINGLE)
1997               fprintf(stderr,"%9.1f-     %5.1f\n",ir->opts.anneal_time[i][j],ir->opts.anneal_temp[i][j]);
1998             else {
1999               fprintf(stderr,"%9.1f      %5.1f\n",ir->opts.anneal_time[i][j],ir->opts.anneal_temp[i][j]);
2000               if(fabs(ir->opts.anneal_temp[i][j]-ir->opts.anneal_temp[i][0])>GMX_REAL_EPS)
2001                 warning_note(wi,"There is a temperature jump when your annealing loops back.\n");
2002             }
2003           }
2004         }
2005       }
2006     }
2007   }
2008
2009   if (ir->ePull != epullNO) {
2010     make_pull_groups(ir->pull,pull_grp,grps,gnames);
2011   }
2012
2013   if (ir->bRot) {
2014     make_rotation_groups(ir->rot,rot_grp,grps,gnames);
2015   }
2016
2017   nacc = str_nelem(acc,MAXPTR,ptr1);
2018   nacg = str_nelem(accgrps,MAXPTR,ptr2);
2019   if (nacg*DIM != nacc)
2020     gmx_fatal(FARGS,"Invalid Acceleration input: %d groups and %d acc. values",
2021                 nacg,nacc);
2022   do_numbering(natoms,groups,nacg,ptr2,grps,gnames,egcACC,
2023                restnm,egrptpALL_GENREST,bVerbose,wi);
2024   nr = groups->grps[egcACC].nr;
2025   snew(ir->opts.acc,nr);
2026   ir->opts.ngacc=nr;
2027
2028   for(i=k=0; (i<nacg); i++)
2029     for(j=0; (j<DIM); j++,k++)
2030       ir->opts.acc[i][j]=strtod(ptr1[k],NULL);
2031   for( ;(i<nr); i++)
2032     for(j=0; (j<DIM); j++)
2033       ir->opts.acc[i][j]=0;
2034
2035   nfrdim  = str_nelem(frdim,MAXPTR,ptr1);
2036   nfreeze = str_nelem(freeze,MAXPTR,ptr2);
2037   if (nfrdim != DIM*nfreeze)
2038     gmx_fatal(FARGS,"Invalid Freezing input: %d groups and %d freeze values",
2039                 nfreeze,nfrdim);
2040   do_numbering(natoms,groups,nfreeze,ptr2,grps,gnames,egcFREEZE,
2041                restnm,egrptpALL_GENREST,bVerbose,wi);
2042   nr = groups->grps[egcFREEZE].nr;
2043   ir->opts.ngfrz=nr;
2044   snew(ir->opts.nFreeze,nr);
2045   for(i=k=0; (i<nfreeze); i++)
2046     for(j=0; (j<DIM); j++,k++) {
2047       ir->opts.nFreeze[i][j]=(gmx_strncasecmp(ptr1[k],"Y",1)==0);
2048       if (!ir->opts.nFreeze[i][j]) {
2049         if (gmx_strncasecmp(ptr1[k],"N",1) != 0) {
2050           sprintf(warnbuf,"Please use Y(ES) or N(O) for freezedim only "
2051                   "(not %s)", ptr1[k]);
2052           warning(wi,warn_buf);
2053         }
2054       }
2055     }
2056   for( ; (i<nr); i++)
2057     for(j=0; (j<DIM); j++)
2058       ir->opts.nFreeze[i][j]=0;
2059
2060   nenergy=str_nelem(energy,MAXPTR,ptr1);
2061   do_numbering(natoms,groups,nenergy,ptr1,grps,gnames,egcENER,
2062                restnm,egrptpALL_GENREST,bVerbose,wi);
2063   add_wall_energrps(groups,ir->nwall,symtab);
2064   ir->opts.ngener = groups->grps[egcENER].nr;
2065   nvcm=str_nelem(vcm,MAXPTR,ptr1);
2066   bRest =
2067     do_numbering(natoms,groups,nvcm,ptr1,grps,gnames,egcVCM,
2068                  restnm,nvcm==0 ? egrptpALL_GENREST : egrptpPART,bVerbose,wi);
2069   if (bRest) {
2070     warning(wi,"Some atoms are not part of any center of mass motion removal group.\n"
2071             "This may lead to artifacts.\n"
2072             "In most cases one should use one group for the whole system.");
2073   }
2074
2075   /* Now we have filled the freeze struct, so we can calculate NRDF */
2076   calc_nrdf(mtop,ir,gnames);
2077
2078   if (v && NULL) {
2079     real fac,ntot=0;
2080
2081     /* Must check per group! */
2082     for(i=0; (i<ir->opts.ngtc); i++)
2083       ntot += ir->opts.nrdf[i];
2084     if (ntot != (DIM*natoms)) {
2085       fac = sqrt(ntot/(DIM*natoms));
2086       if (bVerbose)
2087         fprintf(stderr,"Scaling velocities by a factor of %.3f to account for constraints\n"
2088                 "and removal of center of mass motion\n",fac);
2089       for(i=0; (i<natoms); i++)
2090         svmul(fac,v[i],v[i]);
2091     }
2092   }
2093
2094   nuser=str_nelem(user1,MAXPTR,ptr1);
2095   do_numbering(natoms,groups,nuser,ptr1,grps,gnames,egcUser1,
2096                restnm,egrptpALL_GENREST,bVerbose,wi);
2097   nuser=str_nelem(user2,MAXPTR,ptr1);
2098   do_numbering(natoms,groups,nuser,ptr1,grps,gnames,egcUser2,
2099                restnm,egrptpALL_GENREST,bVerbose,wi);
2100   nuser=str_nelem(xtc_grps,MAXPTR,ptr1);
2101   do_numbering(natoms,groups,nuser,ptr1,grps,gnames,egcXTC,
2102                restnm,egrptpONE,bVerbose,wi);
2103   nofg = str_nelem(orirefitgrp,MAXPTR,ptr1);
2104   do_numbering(natoms,groups,nofg,ptr1,grps,gnames,egcORFIT,
2105                restnm,egrptpALL_GENREST,bVerbose,wi);
2106
2107   /* QMMM input processing */
2108   nQMg          = str_nelem(QMMM,MAXPTR,ptr1);
2109   nQMmethod     = str_nelem(QMmethod,MAXPTR,ptr2);
2110   nQMbasis      = str_nelem(QMbasis,MAXPTR,ptr3);
2111   if((nQMmethod != nQMg)||(nQMbasis != nQMg)){
2112     gmx_fatal(FARGS,"Invalid QMMM input: %d groups %d basissets"
2113               " and %d methods\n",nQMg,nQMbasis,nQMmethod);
2114   }
2115   /* group rest, if any, is always MM! */
2116   do_numbering(natoms,groups,nQMg,ptr1,grps,gnames,egcQMMM,
2117                restnm,egrptpALL_GENREST,bVerbose,wi);
2118   nr = nQMg; /*atoms->grps[egcQMMM].nr;*/
2119   ir->opts.ngQM = nQMg;
2120   snew(ir->opts.QMmethod,nr);
2121   snew(ir->opts.QMbasis,nr);
2122   for(i=0;i<nr;i++){
2123     /* input consists of strings: RHF CASSCF PM3 .. These need to be
2124      * converted to the corresponding enum in names.c
2125      */
2126     ir->opts.QMmethod[i] = search_QMstring(ptr2[i],eQMmethodNR,
2127                                            eQMmethod_names);
2128     ir->opts.QMbasis[i]  = search_QMstring(ptr3[i],eQMbasisNR,
2129                                            eQMbasis_names);
2130
2131   }
2132   nQMmult   = str_nelem(QMmult,MAXPTR,ptr1);
2133   nQMcharge = str_nelem(QMcharge,MAXPTR,ptr2);
2134   nbSH      = str_nelem(bSH,MAXPTR,ptr3);
2135   snew(ir->opts.QMmult,nr);
2136   snew(ir->opts.QMcharge,nr);
2137   snew(ir->opts.bSH,nr);
2138
2139   for(i=0;i<nr;i++){
2140     ir->opts.QMmult[i]   = strtol(ptr1[i],NULL,10);
2141     ir->opts.QMcharge[i] = strtol(ptr2[i],NULL,10);
2142     ir->opts.bSH[i]      = (gmx_strncasecmp(ptr3[i],"Y",1)==0);
2143   }
2144
2145   nCASelec  = str_nelem(CASelectrons,MAXPTR,ptr1);
2146   nCASorb   = str_nelem(CASorbitals,MAXPTR,ptr2);
2147   snew(ir->opts.CASelectrons,nr);
2148   snew(ir->opts.CASorbitals,nr);
2149   for(i=0;i<nr;i++){
2150     ir->opts.CASelectrons[i]= strtol(ptr1[i],NULL,10);
2151     ir->opts.CASorbitals[i] = strtol(ptr2[i],NULL,10);
2152   }
2153   /* special optimization options */
2154
2155   nbOPT = str_nelem(bOPT,MAXPTR,ptr1);
2156   nbTS = str_nelem(bTS,MAXPTR,ptr2);
2157   snew(ir->opts.bOPT,nr);
2158   snew(ir->opts.bTS,nr);
2159   for(i=0;i<nr;i++){
2160     ir->opts.bOPT[i] = (gmx_strncasecmp(ptr1[i],"Y",1)==0);
2161     ir->opts.bTS[i]  = (gmx_strncasecmp(ptr2[i],"Y",1)==0);
2162   }
2163   nSAon     = str_nelem(SAon,MAXPTR,ptr1);
2164   nSAoff    = str_nelem(SAoff,MAXPTR,ptr2);
2165   nSAsteps  = str_nelem(SAsteps,MAXPTR,ptr3);
2166   snew(ir->opts.SAon,nr);
2167   snew(ir->opts.SAoff,nr);
2168   snew(ir->opts.SAsteps,nr);
2169
2170   for(i=0;i<nr;i++){
2171     ir->opts.SAon[i]    = strtod(ptr1[i],NULL);
2172     ir->opts.SAoff[i]   = strtod(ptr2[i],NULL);
2173     ir->opts.SAsteps[i] = strtol(ptr3[i],NULL,10);
2174   }
2175   /* end of QMMM input */
2176
2177   if (bVerbose)
2178     for(i=0; (i<egcNR); i++) {
2179       fprintf(stderr,"%-16s has %d element(s):",gtypes[i],groups->grps[i].nr);
2180       for(j=0; (j<groups->grps[i].nr); j++)
2181         fprintf(stderr," %s",*(groups->grpname[groups->grps[i].nm_ind[j]]));
2182       fprintf(stderr,"\n");
2183     }
2184
2185   nr = groups->grps[egcENER].nr;
2186   snew(ir->opts.egp_flags,nr*nr);
2187
2188   bExcl = do_egp_flag(ir,groups,"energygrp-excl",egpexcl,EGP_EXCL);
2189   if (bExcl && EEL_FULL(ir->coulombtype))
2190     warning(wi,"Can not exclude the lattice Coulomb energy between energy groups");
2191
2192   bTable = do_egp_flag(ir,groups,"energygrp-table",egptable,EGP_TABLE);
2193   if (bTable && !(ir->vdwtype == evdwUSER) &&
2194       !(ir->coulombtype == eelUSER) && !(ir->coulombtype == eelPMEUSER) &&
2195       !(ir->coulombtype == eelPMEUSERSWITCH))
2196     gmx_fatal(FARGS,"Can only have energy group pair tables in combination with user tables for VdW and/or Coulomb");
2197
2198   decode_cos(efield_x,&(ir->ex[XX]),FALSE);
2199   decode_cos(efield_xt,&(ir->et[XX]),TRUE);
2200   decode_cos(efield_y,&(ir->ex[YY]),FALSE);
2201   decode_cos(efield_yt,&(ir->et[YY]),TRUE);
2202   decode_cos(efield_z,&(ir->ex[ZZ]),FALSE);
2203   decode_cos(efield_zt,&(ir->et[ZZ]),TRUE);
2204
2205   for(i=0; (i<grps->nr); i++)
2206     sfree(gnames[i]);
2207   sfree(gnames);
2208   done_blocka(grps);
2209   sfree(grps);
2210
2211 }
2212
2213
2214
2215 static void check_disre(gmx_mtop_t *mtop)
2216 {
2217   gmx_ffparams_t *ffparams;
2218   t_functype *functype;
2219   t_iparams  *ip;
2220   int i,ndouble,ftype;
2221   int label,old_label;
2222
2223   if (gmx_mtop_ftype_count(mtop,F_DISRES) > 0) {
2224     ffparams  = &mtop->ffparams;
2225     functype  = ffparams->functype;
2226     ip        = ffparams->iparams;
2227     ndouble   = 0;
2228     old_label = -1;
2229     for(i=0; i<ffparams->ntypes; i++) {
2230       ftype = functype[i];
2231       if (ftype == F_DISRES) {
2232         label = ip[i].disres.label;
2233         if (label == old_label) {
2234           fprintf(stderr,"Distance restraint index %d occurs twice\n",label);
2235           ndouble++;
2236         }
2237         old_label = label;
2238       }
2239     }
2240     if (ndouble>0)
2241       gmx_fatal(FARGS,"Found %d double distance restraint indices,\n"
2242                 "probably the parameters for multiple pairs in one restraint "
2243                 "are not identical\n",ndouble);
2244   }
2245 }
2246
2247 static gmx_bool absolute_reference(t_inputrec *ir,gmx_mtop_t *sys,
2248                                    gmx_bool posres_only,
2249                                    ivec AbsRef)
2250 {
2251     int d,g,i;
2252     gmx_mtop_ilistloop_t iloop;
2253     t_ilist *ilist;
2254     int nmol;
2255     t_iparams *pr;
2256
2257     clear_ivec(AbsRef);
2258
2259     if (!posres_only)
2260     {
2261         /* Check the COM */
2262         for(d=0; d<DIM; d++)
2263         {
2264             AbsRef[d] = (d < ndof_com(ir) ? 0 : 1);
2265         }
2266         /* Check for freeze groups */
2267         for(g=0; g<ir->opts.ngfrz; g++)
2268         {
2269             for(d=0; d<DIM; d++)
2270             {
2271                 if (ir->opts.nFreeze[g][d] != 0)
2272                 {
2273                     AbsRef[d] = 1;
2274                 }
2275             }
2276         }
2277     }
2278
2279     /* Check for position restraints */
2280     iloop = gmx_mtop_ilistloop_init(sys);
2281     while (gmx_mtop_ilistloop_next(iloop,&ilist,&nmol))
2282     {
2283         if (nmol > 0 &&
2284             (AbsRef[XX] == 0 || AbsRef[YY] == 0 || AbsRef[ZZ] == 0))
2285         {
2286             for(i=0; i<ilist[F_POSRES].nr; i+=2)
2287             {
2288                 pr = &sys->ffparams.iparams[ilist[F_POSRES].iatoms[i]];
2289                 for(d=0; d<DIM; d++)
2290                 {
2291                     if (pr->posres.fcA[d] != 0)
2292                     {
2293                         AbsRef[d] = 1;
2294                     }
2295                 }
2296             }
2297         }
2298     }
2299
2300     return (AbsRef[XX] != 0 && AbsRef[YY] != 0 && AbsRef[ZZ] != 0);
2301 }
2302
2303 void triple_check(const char *mdparin,t_inputrec *ir,gmx_mtop_t *sys,
2304                   warninp_t wi)
2305 {
2306   char err_buf[256];
2307   int  i,m,g,nmol,npct;
2308   gmx_bool bCharge,bAcc;
2309   real gdt_max,*mgrp,mt;
2310   rvec acc;
2311   gmx_mtop_atomloop_block_t aloopb;
2312   gmx_mtop_atomloop_all_t aloop;
2313   t_atom *atom;
2314   ivec AbsRef;
2315   char warn_buf[STRLEN];
2316
2317   set_warning_line(wi,mdparin,-1);
2318
2319   if (EI_DYNAMICS(ir->eI) && !EI_SD(ir->eI) && ir->eI != eiBD &&
2320       ir->comm_mode == ecmNO &&
2321       !(absolute_reference(ir,sys,FALSE,AbsRef) || ir->nsteps <= 10)) {
2322     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");
2323   }
2324
2325     /* Check for pressure coupling with absolute position restraints */
2326     if (ir->epc != epcNO && ir->refcoord_scaling == erscNO)
2327     {
2328         absolute_reference(ir,sys,TRUE,AbsRef);
2329         {
2330             for(m=0; m<DIM; m++)
2331             {
2332                 if (AbsRef[m] && norm2(ir->compress[m]) > 0)
2333                 {
2334                     warning(wi,"You are using pressure coupling with absolute position restraints, this will give artifacts. Use the refcoord_scaling option.");
2335                     break;
2336                 }
2337             }
2338         }
2339     }
2340
2341   bCharge = FALSE;
2342   aloopb = gmx_mtop_atomloop_block_init(sys);
2343   while (gmx_mtop_atomloop_block_next(aloopb,&atom,&nmol)) {
2344     if (atom->q != 0 || atom->qB != 0) {
2345       bCharge = TRUE;
2346     }
2347   }
2348
2349   if (!bCharge) {
2350     if (EEL_FULL(ir->coulombtype)) {
2351       sprintf(err_buf,
2352               "You are using full electrostatics treatment %s for a system without charges.\n"
2353               "This costs a lot of performance for just processing zeros, consider using %s instead.\n",
2354               EELTYPE(ir->coulombtype),EELTYPE(eelCUT));
2355       warning(wi,err_buf);
2356     }
2357   } else {
2358     if (ir->coulombtype == eelCUT && ir->rcoulomb > 0 && !ir->implicit_solvent) {
2359       sprintf(err_buf,
2360               "You are using a plain Coulomb cut-off, which might produce artifacts.\n"
2361               "You might want to consider using %s electrostatics.\n",
2362               EELTYPE(eelPME));
2363       warning_note(wi,err_buf);
2364     }
2365   }
2366
2367   /* Generalized reaction field */
2368   if (ir->opts.ngtc == 0) {
2369     sprintf(err_buf,"No temperature coupling while using coulombtype %s",
2370             eel_names[eelGRF]);
2371     CHECK(ir->coulombtype == eelGRF);
2372   }
2373   else {
2374     sprintf(err_buf,"When using coulombtype = %s"
2375             " ref_t for temperature coupling should be > 0",
2376             eel_names[eelGRF]);
2377     CHECK((ir->coulombtype == eelGRF) && (ir->opts.ref_t[0] <= 0));
2378   }
2379
2380   if (ir->eI == eiSD1) {
2381     gdt_max = 0;
2382     for(i=0; (i<ir->opts.ngtc); i++)
2383       gdt_max = max(gdt_max,ir->delta_t/ir->opts.tau_t[i]);
2384     if (0.5*gdt_max > 0.0015) {
2385       sprintf(warn_buf,"The relative error with integrator %s is 0.5*delta-t/tau-t = %g, you might want to switch to integrator %s\n",
2386               ei_names[ir->eI],0.5*gdt_max,ei_names[eiSD2]);
2387       warning_note(wi,warn_buf);
2388     }
2389   }
2390
2391   bAcc = FALSE;
2392   for(i=0; (i<sys->groups.grps[egcACC].nr); i++) {
2393     for(m=0; (m<DIM); m++) {
2394       if (fabs(ir->opts.acc[i][m]) > 1e-6) {
2395         bAcc = TRUE;
2396       }
2397     }
2398   }
2399   if (bAcc) {
2400     clear_rvec(acc);
2401     snew(mgrp,sys->groups.grps[egcACC].nr);
2402     aloop = gmx_mtop_atomloop_all_init(sys);
2403     while (gmx_mtop_atomloop_all_next(aloop,&i,&atom)) {
2404       mgrp[ggrpnr(&sys->groups,egcACC,i)] += atom->m;
2405     }
2406     mt = 0.0;
2407     for(i=0; (i<sys->groups.grps[egcACC].nr); i++) {
2408       for(m=0; (m<DIM); m++)
2409         acc[m] += ir->opts.acc[i][m]*mgrp[i];
2410       mt += mgrp[i];
2411     }
2412     for(m=0; (m<DIM); m++) {
2413       if (fabs(acc[m]) > 1e-6) {
2414         const char *dim[DIM] = { "X", "Y", "Z" };
2415         fprintf(stderr,
2416                 "Net Acceleration in %s direction, will %s be corrected\n",
2417                 dim[m],ir->nstcomm != 0 ? "" : "not");
2418         if (ir->nstcomm != 0 && m < ndof_com(ir)) {
2419           acc[m] /= mt;
2420           for (i=0; (i<sys->groups.grps[egcACC].nr); i++)
2421             ir->opts.acc[i][m] -= acc[m];
2422         }
2423       }
2424     }
2425     sfree(mgrp);
2426   }
2427
2428   if (ir->efep != efepNO && ir->sc_alpha != 0 &&
2429       !gmx_within_tol(sys->ffparams.reppow,12.0,10*GMX_DOUBLE_EPS)) {
2430     gmx_fatal(FARGS,"Soft-core interactions are only supported with VdW repulsion power 12");
2431   }
2432
2433   if (ir->ePull != epullNO) {
2434     if (ir->pull->grp[0].nat == 0) {
2435         absolute_reference(ir,sys,FALSE,AbsRef);
2436       for(m=0; m<DIM; m++) {
2437         if (ir->pull->dim[m] && !AbsRef[m]) {
2438           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.");
2439           break;
2440         }
2441       }
2442     }
2443
2444     if (ir->pull->eGeom == epullgDIRPBC) {
2445       for(i=0; i<3; i++) {
2446         for(m=0; m<=i; m++) {
2447           if ((ir->epc != epcNO && ir->compress[i][m] != 0) ||
2448               ir->deform[i][m] != 0) {
2449             for(g=1; g<ir->pull->ngrp; g++) {
2450               if (ir->pull->grp[g].vec[m] != 0) {
2451                 gmx_fatal(FARGS,"Can not have dynamic box while using pull geometry '%s' (dim %c)",EPULLGEOM(ir->pull->eGeom),'x'+m);
2452               }
2453             }
2454           }
2455         }
2456       }
2457     }
2458   }
2459
2460   check_disre(sys);
2461 }
2462
2463 void double_check(t_inputrec *ir,matrix box,gmx_bool bConstr,warninp_t wi)
2464 {
2465   real min_size;
2466   gmx_bool bTWIN;
2467   char warn_buf[STRLEN];
2468   const char *ptr;
2469
2470   ptr = check_box(ir->ePBC,box);
2471   if (ptr) {
2472       warning_error(wi,ptr);
2473   }
2474
2475   if (bConstr && ir->eConstrAlg == econtSHAKE) {
2476     if (ir->shake_tol <= 0.0) {
2477       sprintf(warn_buf,"ERROR: shake-tol must be > 0 instead of %g\n",
2478               ir->shake_tol);
2479       warning_error(wi,warn_buf);
2480     }
2481
2482     if (IR_TWINRANGE(*ir) && ir->nstlist > 1) {
2483       sprintf(warn_buf,"With twin-range cut-off's and SHAKE the virial and the pressure are incorrect.");
2484       if (ir->epc == epcNO) {
2485         warning(wi,warn_buf);
2486       } else {
2487           warning_error(wi,warn_buf);
2488       }
2489     }
2490   }
2491
2492   if( (ir->eConstrAlg == econtLINCS) && bConstr) {
2493     /* If we have Lincs constraints: */
2494     if(ir->eI==eiMD && ir->etc==etcNO &&
2495        ir->eConstrAlg==econtLINCS && ir->nLincsIter==1) {
2496       sprintf(warn_buf,"For energy conservation with LINCS, lincs_iter should be 2 or larger.\n");
2497       warning_note(wi,warn_buf);
2498     }
2499
2500     if ((ir->eI == eiCG || ir->eI == eiLBFGS) && (ir->nProjOrder<8)) {
2501       sprintf(warn_buf,"For accurate %s with LINCS constraints, lincs-order should be 8 or more.",ei_names[ir->eI]);
2502       warning_note(wi,warn_buf);
2503     }
2504     if (ir->epc==epcMTTK) {
2505         warning_error(wi,"MTTK not compatible with lincs -- use shake instead.");
2506     }
2507   }
2508
2509   if (ir->LincsWarnAngle > 90.0) {
2510     sprintf(warn_buf,"lincs-warnangle can not be larger than 90 degrees, setting it to 90.\n");
2511     warning(wi,warn_buf);
2512     ir->LincsWarnAngle = 90.0;
2513   }
2514
2515   if (ir->ePBC != epbcNONE) {
2516     if (ir->nstlist == 0) {
2517       warning(wi,"With nstlist=0 atoms are only put into the box at step 0, therefore drifting atoms might cause the simulation to crash.");
2518     }
2519     bTWIN = (ir->rlistlong > ir->rlist);
2520     if (ir->ns_type == ensGRID) {
2521       if (sqr(ir->rlistlong) >= max_cutoff2(ir->ePBC,box)) {
2522           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",
2523                 bTWIN ? (ir->rcoulomb==ir->rlistlong ? "rcoulomb" : "rvdw"):"rlist");
2524           warning_error(wi,warn_buf);
2525       }
2526     } else {
2527       min_size = min(box[XX][XX],min(box[YY][YY],box[ZZ][ZZ]));
2528       if (2*ir->rlistlong >= min_size) {
2529           sprintf(warn_buf,"ERROR: One of the box lengths is smaller than twice the cut-off length. Increase the box size or decrease rlist.");
2530           warning_error(wi,warn_buf);
2531         if (TRICLINIC(box))
2532           fprintf(stderr,"Grid search might allow larger cut-off's than simple search with triclinic boxes.");
2533       }
2534     }
2535   }
2536 }
2537
2538 void check_chargegroup_radii(const gmx_mtop_t *mtop,const t_inputrec *ir,
2539                              rvec *x,
2540                              warninp_t wi)
2541 {
2542     real rvdw1,rvdw2,rcoul1,rcoul2;
2543     char warn_buf[STRLEN];
2544
2545     calc_chargegroup_radii(mtop,x,&rvdw1,&rvdw2,&rcoul1,&rcoul2);
2546
2547     if (rvdw1 > 0)
2548     {
2549         printf("Largest charge group radii for Van der Waals: %5.3f, %5.3f nm\n",
2550                rvdw1,rvdw2);
2551     }
2552     if (rcoul1 > 0)
2553     {
2554         printf("Largest charge group radii for Coulomb:       %5.3f, %5.3f nm\n",
2555                rcoul1,rcoul2);
2556     }
2557
2558     if (ir->rlist > 0)
2559     {
2560         if (rvdw1  + rvdw2  > ir->rlist ||
2561             rcoul1 + rcoul2 > ir->rlist)
2562         {
2563             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);
2564             warning(wi,warn_buf);
2565         }
2566         else
2567         {
2568             /* Here we do not use the zero at cut-off macro,
2569              * since user defined interactions might purposely
2570              * not be zero at the cut-off.
2571              */
2572             if (EVDW_IS_ZERO_AT_CUTOFF(ir->vdwtype) &&
2573                 rvdw1 + rvdw2 > ir->rlist - ir->rvdw)
2574             {
2575                 sprintf(warn_buf,"The sum of the two largest charge group radii (%f) is larger than rlist (%f) - rvdw (%f)\n",
2576                         rvdw1+rvdw2,
2577                         ir->rlist,ir->rvdw);
2578                 if (ir_NVE(ir))
2579                 {
2580                     warning(wi,warn_buf);
2581                 }
2582                 else
2583                 {
2584                     warning_note(wi,warn_buf);
2585                 }
2586             }
2587             if (EEL_IS_ZERO_AT_CUTOFF(ir->coulombtype) &&
2588                 rcoul1 + rcoul2 > ir->rlistlong - ir->rcoulomb)
2589             {
2590                 sprintf(warn_buf,"The sum of the two largest charge group radii (%f) is larger than %s (%f) - rcoulomb (%f)\n",
2591                         rcoul1+rcoul2,
2592                         ir->rlistlong > ir->rlist ? "rlistlong" : "rlist",
2593                         ir->rlistlong,ir->rcoulomb);
2594                 if (ir_NVE(ir))
2595                 {
2596                     warning(wi,warn_buf);
2597                 }
2598                 else
2599                 {
2600                     warning_note(wi,warn_buf);
2601                 }
2602             }
2603         }
2604     }
2605 }