src/gromacs/gmxpreprocess/readir.c

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