src/kernel/readir.c

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