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