src/gromacs/gmxpreprocess/readir.c

   1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
   2  *
   3  *
   4  *                This source code is part of
   5  *
   6  *                 G   R   O   M   A   C   S
   7  *
   8  *          GROningen MAchine for Chemical Simulations
   9  *
  10  *                        VERSION 3.2.0
  11  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  12  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  13  * Copyright (c) 2001-2004, The GROMACS development team,
  14  * check out http://www.gromacs.org for more information.
  15
  16  * This program is free software; you can redistribute it and/or
  17  * modify it under the terms of the GNU General Public License
  18  * as published by the Free Software Foundation; either version 2
  19  * of the License, or (at your option) any later version.
  20  *
  21  * If you want to redistribute modifications, please consider that
  22  * scientific software is very special. Version control is crucial -
  23  * bugs must be traceable. We will be happy to consider code for
  24  * inclusion in the official distribution, but derived work must not
  25  * be called official GROMACS. Details are found in the README & COPYING
  26  * files - if they are missing, get the official version at www.gromacs.org.
  27  *
  28  * To help us fund GROMACS development, we humbly ask that you cite
  29  * the papers on the package - you can find them in the top README file.
  30  *
  31  * For more info, check our website at http://www.gromacs.org
  32  *
  33  * And Hey:
  34  * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
  35  */
  36 #ifdef HAVE_CONFIG_H
  37 #include <config.h>
  38 #endif
  39
  40 #include <ctype.h>
  41 #include <stdlib.h>
  42 #include <limits.h>
  43 #include "sysstuff.h"
  44 #include "smalloc.h"
  45 #include "typedefs.h"
  46 #include "physics.h"
  47 #include "names.h"
  48 #include "gmx_fatal.h"
  49 #include "macros.h"
  50 #include "index.h"
  51 #include "symtab.h"
  52 #include "string2.h"
  53 #include "readinp.h"
  54 #include "warninp.h"
  55 #include "readir.h"
  56 #include "toputil.h"
  57 #include "index.h"
  58 #include "network.h"
  59 #include "vec.h"
  60 #include "pbc.h"
  61 #include "mtop_util.h"
  62 #include "chargegroup.h"
  63 #include "inputrec.h"
  64
  65 #define MAXPTR 254
  66 #define NOGID  255
  67 #define MAXLAMBDAS 1024
  68
  69 /* Resource parameters
  70  * Do not change any of these until you read the instruction
  71  * in readinp.h. Some cpp's do not take spaces after the backslash
  72  * (like the c-shell), which will give you a very weird compiler
  73  * message.
  74  */
  75
  76 static char tcgrps[STRLEN], tau_t[STRLEN], ref_t[STRLEN],
  77             acc[STRLEN], accgrps[STRLEN], freeze[STRLEN], frdim[STRLEN],
  78             energy[STRLEN], user1[STRLEN], user2[STRLEN], vcm[STRLEN], xtc_grps[STRLEN],
  79             couple_moltype[STRLEN], orirefitgrp[STRLEN], egptable[STRLEN], egpexcl[STRLEN],
  80             wall_atomtype[STRLEN], wall_density[STRLEN], deform[STRLEN], QMMM[STRLEN];
  81 static char   fep_lambda[efptNR][STRLEN];
  82 static char   lambda_weights[STRLEN];
  83 static char **pull_grp;
  84 static char **rot_grp;
  85 static char   anneal[STRLEN], anneal_npoints[STRLEN],
  86               anneal_time[STRLEN], anneal_temp[STRLEN];
  87 static char   QMmethod[STRLEN], QMbasis[STRLEN], QMcharge[STRLEN], QMmult[STRLEN],
  88               bSH[STRLEN], CASorbitals[STRLEN], CASelectrons[STRLEN], SAon[STRLEN],
  89               SAoff[STRLEN], SAsteps[STRLEN], bTS[STRLEN], bOPT[STRLEN];
  90 static char efield_x[STRLEN], efield_xt[STRLEN], efield_y[STRLEN],
  91             efield_yt[STRLEN], efield_z[STRLEN], efield_zt[STRLEN];
  92
  93 enum {
  94     egrptpALL,         /* All particles have to be a member of a group.     */
  95     egrptpALL_GENREST, /* A rest group with name is generated for particles *
  96                         * that are not part of any group.                   */
  97     egrptpPART,        /* As egrptpALL_GENREST, but no name is generated    *
  98                         * for the rest group.                               */
  99     egrptpONE          /* Merge all selected groups into one group,         *
 100                         * make a rest group for the remaining particles.    */
 101 };
 102
 103 static const char *constraints[eshNR+1]    = {
 104     "none", "h-bonds", "all-bonds", "h-angles", "all-angles", NULL
 105 };
 106
 107 static const char *couple_lam[ecouplamNR+1]    = {
 108     "vdw-q", "vdw", "q", "none", NULL
 109 };
 110
 111 void init_ir(t_inputrec *ir, t_gromppopts *opts)
 112 {
 113     snew(opts->include, STRLEN);
 114     snew(opts->define, STRLEN);
 115     snew(ir->fepvals, 1);
 116     snew(ir->expandedvals, 1);
 117     snew(ir->simtempvals, 1);
 118 }
 119
 120 static void GetSimTemps(int ntemps, t_simtemp *simtemp, double *temperature_lambdas)
 121 {
 122
 123     int i;
 124
 125     for (i = 0; i < ntemps; i++)
 126     {
 127         /* simple linear scaling -- allows more control */
 128         if (simtemp->eSimTempScale == esimtempLINEAR)
 129         {
 130             simtemp->temperatures[i] = simtemp->simtemp_low + (simtemp->simtemp_high-simtemp->simtemp_low)*temperature_lambdas[i];
 131         }
 132         else if (simtemp->eSimTempScale == esimtempGEOMETRIC)  /* should give roughly equal acceptance for constant heat capacity . . . */
 133         {
 134             simtemp->temperatures[i] = simtemp->simtemp_low * pow(simtemp->simtemp_high/simtemp->simtemp_low, (1.0*i)/(ntemps-1));
 135         }
 136         else if (simtemp->eSimTempScale == esimtempEXPONENTIAL)
 137         {
 138             simtemp->temperatures[i] = simtemp->simtemp_low + (simtemp->simtemp_high-simtemp->simtemp_low)*((exp(temperature_lambdas[i])-1)/(exp(1.0)-1));
 139         }
 140         else
 141         {
 142             char errorstr[128];
 143             sprintf(errorstr, "eSimTempScale=%d not defined", simtemp->eSimTempScale);
 144             gmx_fatal(FARGS, errorstr);
 145         }
 146     }
 147 }
 148
 149
 150
 151 static void _low_check(gmx_bool b, char *s, warninp_t wi)
 152 {
 153     if (b)
 154     {
 155         warning_error(wi, s);
 156     }
 157 }
 158
 159 static void check_nst(const char *desc_nst, int nst,
 160                       const char *desc_p, int *p,
 161                       warninp_t wi)
 162 {
 163     char buf[STRLEN];
 164
 165     if (*p > 0 && *p % nst != 0)
 166     {
 167         /* Round up to the next multiple of nst */
 168         *p = ((*p)/nst + 1)*nst;
 169         sprintf(buf, "%s should be a multiple of %s, changing %s to %d\n",
 170                 desc_p, desc_nst, desc_p, *p);
 171         warning(wi, buf);
 172     }
 173 }
 174
 175 static gmx_bool ir_NVE(const t_inputrec *ir)
 176 {
 177     return ((ir->eI == eiMD || EI_VV(ir->eI)) && ir->etc == etcNO);
 178 }
 179
 180 static int lcd(int n1, int n2)
 181 {
 182     int d, i;
 183
 184     d = 1;
 185     for (i = 2; (i <= n1 && i <= n2); i++)
 186     {
 187         if (n1 % i == 0 && n2 % i == 0)
 188         {
 189             d = i;
 190         }
 191     }
 192
 193     return d;
 194 }
 195
 196 static void process_interaction_modifier(const t_inputrec *ir, int *eintmod)
 197 {
 198     if (*eintmod == eintmodPOTSHIFT_VERLET)
 199     {
 200         if (ir->cutoff_scheme == ecutsVERLET)
 201         {
 202             *eintmod = eintmodPOTSHIFT;
 203         }
 204         else
 205         {
 206             *eintmod = eintmodNONE;
 207         }
 208     }
 209 }
 210
 211 void check_ir(const char *mdparin, t_inputrec *ir, t_gromppopts *opts,
 212               warninp_t wi)
 213 /* Check internal consistency */
 214 {
 215     /* Strange macro: first one fills the err_buf, and then one can check
 216      * the condition, which will print the message and increase the error
 217      * counter.
 218      */
 219 #define CHECK(b) _low_check(b, err_buf, wi)
 220     char        err_buf[256], warn_buf[STRLEN];
 221     int         i, j;
 222     int         ns_type  = 0;
 223     real        dt_coupl = 0;
 224     real        dt_pcoupl;
 225     int         nstcmin;
 226     t_lambda   *fep    = ir->fepvals;
 227     t_expanded *expand = ir->expandedvals;
 228
 229     set_warning_line(wi, mdparin, -1);
 230
 231     /* BASIC CUT-OFF STUFF */
 232     if (ir->rcoulomb < 0)
 233     {
 234         warning_error(wi, "rcoulomb should be >= 0");
 235     }
 236     if (ir->rvdw < 0)
 237     {
 238         warning_error(wi, "rvdw should be >= 0");
 239     }
 240     if (ir->rlist < 0 &&
 241         !(ir->cutoff_scheme == ecutsVERLET && ir->verletbuf_drift > 0))
 242     {
 243         warning_error(wi, "rlist should be >= 0");
 244     }
 245
 246     process_interaction_modifier(ir, &ir->coulomb_modifier);
 247     process_interaction_modifier(ir, &ir->vdw_modifier);
 248
 249     if (ir->cutoff_scheme == ecutsGROUP)
 250     {
 251         /* BASIC CUT-OFF STUFF */
 252         if (ir->rlist == 0 ||
 253             !((EEL_MIGHT_BE_ZERO_AT_CUTOFF(ir->coulombtype) && ir->rcoulomb > ir->rlist) ||
 254               (EVDW_MIGHT_BE_ZERO_AT_CUTOFF(ir->vdwtype)    && ir->rvdw     > ir->rlist)))
 255         {
 256             /* No switched potential and/or no twin-range:
 257              * we can set the long-range cut-off to the maximum of the other cut-offs.
 258              */
 259             ir->rlistlong = max_cutoff(ir->rlist, max_cutoff(ir->rvdw, ir->rcoulomb));
 260         }
 261         else if (ir->rlistlong < 0)
 262         {
 263             ir->rlistlong = max_cutoff(ir->rlist, max_cutoff(ir->rvdw, ir->rcoulomb));
 264             sprintf(warn_buf, "rlistlong was not set, setting it to %g (no buffer)",
 265                     ir->rlistlong);
 266             warning(wi, warn_buf);
 267         }
 268         if (ir->rlistlong == 0 && ir->ePBC != epbcNONE)
 269         {
 270             warning_error(wi, "Can not have an infinite cut-off with PBC");
 271         }
 272         if (ir->rlistlong > 0 && (ir->rlist == 0 || ir->rlistlong < ir->rlist))
 273         {
 274             warning_error(wi, "rlistlong can not be shorter than rlist");
 275         }
 276         if (IR_TWINRANGE(*ir) && ir->nstlist <= 0)
 277         {
 278             warning_error(wi, "Can not have nstlist<=0 with twin-range interactions");
 279         }
 280     }
 281
 282     if (ir->rlistlong == ir->rlist)
 283     {
 284         ir->nstcalclr = 0;
 285     }
 286     else if (ir->rlistlong > ir->rlist && ir->nstcalclr == 0)
 287     {
 288         warning_error(wi, "With different cutoffs for electrostatics and VdW, nstcalclr must be -1 or a positive number");
 289     }
 290
 291     if (ir->cutoff_scheme == ecutsVERLET)
 292     {
 293         real rc_max;
 294
 295         /* Normal Verlet type neighbor-list, currently only limited feature support */
 296         if (inputrec2nboundeddim(ir) < 3)
 297         {
 298             warning_error(wi, "With Verlet lists only full pbc or pbc=xy with walls is supported");
 299         }
 300         if (ir->rcoulomb != ir->rvdw)
 301         {
 302             warning_error(wi, "With Verlet lists rcoulomb!=rvdw is not supported");
 303         }
 304         if (ir->vdwtype != evdwCUT)
 305         {
 306             warning_error(wi, "With Verlet lists only cut-off LJ interactions are supported");
 307         }
 308         if (!(ir->coulombtype == eelCUT ||
 309               (EEL_RF(ir->coulombtype) && ir->coulombtype != eelRF_NEC) ||
 310               EEL_PME(ir->coulombtype) || ir->coulombtype == eelEWALD))
 311         {
 312             warning_error(wi, "With Verlet lists only cut-off, reaction-field, PME and Ewald electrostatics are supported");
 313         }
 314
 315         if (ir->nstlist <= 0)
 316         {
 317             warning_error(wi, "With Verlet lists nstlist should be larger than 0");
 318         }
 319
 320         if (ir->nstlist < 10)
 321         {
 322             warning_note(wi, "With Verlet lists the optimal nstlist is >= 10, with GPUs >= 20. Note that with the Verlet scheme, nstlist has no effect on the accuracy of your simulation.");
 323         }
 324
 325         rc_max = max(ir->rvdw, ir->rcoulomb);
 326
 327         if (ir->verletbuf_drift <= 0)
 328         {
 329             if (ir->verletbuf_drift == 0)
 330             {
 331                 warning_error(wi, "Can not have an energy drift of exactly 0");
 332             }
 333
 334             if (ir->rlist < rc_max)
 335             {
 336                 warning_error(wi, "With verlet lists rlist can not be smaller than rvdw or rcoulomb");
 337             }
 338
 339             if (ir->rlist == rc_max && ir->nstlist > 1)
 340             {
 341                 warning_note(wi, "rlist is equal to rvdw and/or rcoulomb: there is no explicit Verlet buffer. The cluster pair list does have a buffering effect, but choosing a larger rlist might be necessary for good energy conservation.");
 342             }
 343         }
 344         else
 345         {
 346             if (ir->rlist > rc_max)
 347             {
 348                 warning_note(wi, "You have set rlist larger than the interaction cut-off, but you also have verlet-buffer-drift > 0. Will set rlist using verlet-buffer-drift.");
 349             }
 350
 351             if (ir->nstlist == 1)
 352             {
 353                 /* No buffer required */
 354                 ir->rlist = rc_max;
 355             }
 356             else
 357             {
 358                 if (EI_DYNAMICS(ir->eI))
 359                 {
 360                     if (EI_MD(ir->eI) && ir->etc == etcNO)
 361                     {
 362                         warning_error(wi, "Temperature coupling is required for calculating rlist using the energy drift with verlet-buffer-drift > 0. Either use temperature coupling or set rlist yourself together with verlet-buffer-drift = -1.");
 363                     }
 364
 365                     if (inputrec2nboundeddim(ir) < 3)
 366                     {
 367                         warning_error(wi, "The box volume is required for calculating rlist from the energy drift with verlet-buffer-drift > 0. You are using at least one unbounded dimension, so no volume can be computed. Either use a finite box, or set rlist yourself together with verlet-buffer-drift = -1.");
 368                     }
 369                     /* Set rlist temporarily so we can continue processing */
 370                     ir->rlist = rc_max;
 371                 }
 372                 else
 373                 {
 374                     /* Set the buffer to 5% of the cut-off */
 375                     ir->rlist = 1.05*rc_max;
 376                 }
 377             }
 378         }
 379
 380         /* No twin-range calculations with Verlet lists */
 381         ir->rlistlong = ir->rlist;
 382     }
 383
 384     if (ir->nstcalclr == -1)
 385     {
 386         /* if rlist=rlistlong, this will later be changed to nstcalclr=0 */
 387         ir->nstcalclr = ir->nstlist;
 388     }
 389     else if (ir->nstcalclr > 0)
 390     {
 391         if (ir->nstlist > 0 && (ir->nstlist % ir->nstcalclr != 0))
 392         {
 393             warning_error(wi, "nstlist must be evenly divisible by nstcalclr. Use nstcalclr = -1 to automatically follow nstlist");
 394         }
 395     }
 396     else if (ir->nstcalclr < -1)
 397     {
 398         warning_error(wi, "nstcalclr must be a positive number (divisor of nstcalclr), or -1 to follow nstlist.");
 399     }
 400
 401     if (EEL_PME(ir->coulombtype) && ir->rcoulomb > ir->rvdw && ir->nstcalclr > 1)
 402     {
 403         warning_error(wi, "When used with PME, the long-range component of twin-range interactions must be updated every step (nstcalclr)");
 404     }
 405
 406     /* GENERAL INTEGRATOR STUFF */
 407     if (!(ir->eI == eiMD || EI_VV(ir->eI)))
 408     {
 409         ir->etc = etcNO;
 410     }
 411     if (ir->eI == eiVVAK)
 412     {
 413         sprintf(warn_buf, "Integrator method %s is implemented primarily for validation purposes; for molecular dynamics, you should probably be using %s or %s", ei_names[eiVVAK], ei_names[eiMD], ei_names[eiVV]);
 414         warning_note(wi, warn_buf);
 415     }
 416     if (!EI_DYNAMICS(ir->eI))
 417     {
 418         ir->epc = epcNO;
 419     }
 420     if (EI_DYNAMICS(ir->eI))
 421     {
 422         if (ir->nstcalcenergy < 0)
 423         {
 424             ir->nstcalcenergy = ir_optimal_nstcalcenergy(ir);
 425             if (ir->nstenergy != 0 && ir->nstenergy < ir->nstcalcenergy)
 426             {
 427                 /* nstcalcenergy larger than nstener does not make sense.
 428                  * We ideally want nstcalcenergy=nstener.
 429                  */
 430                 if (ir->nstlist > 0)
 431                 {
 432                     ir->nstcalcenergy = lcd(ir->nstenergy, ir->nstlist);
 433                 }
 434                 else
 435                 {
 436                     ir->nstcalcenergy = ir->nstenergy;
 437                 }
 438             }
 439         }
 440         else if ( (ir->nstenergy > 0 && ir->nstcalcenergy > ir->nstenergy) ||
 441                   (ir->efep != efepNO && ir->fepvals->nstdhdl > 0 &&
 442                    (ir->nstcalcenergy > ir->fepvals->nstdhdl) ) )
 443
 444         {
 445             const char *nsten    = "nstenergy";
 446             const char *nstdh    = "nstdhdl";
 447             const char *min_name = nsten;
 448             int         min_nst  = ir->nstenergy;
 449
 450             /* find the smallest of ( nstenergy, nstdhdl ) */
 451             if (ir->efep != efepNO && ir->fepvals->nstdhdl > 0 &&
 452                 (ir->fepvals->nstdhdl < ir->nstenergy) )
 453             {
 454                 min_nst  = ir->fepvals->nstdhdl;
 455                 min_name = nstdh;
 456             }
 457             /* If the user sets nstenergy small, we should respect that */
 458             sprintf(warn_buf,
 459                     "Setting nstcalcenergy (%d) equal to %s (%d)",
 460                     ir->nstcalcenergy, min_name, min_nst);
 461             warning_note(wi, warn_buf);
 462             ir->nstcalcenergy = min_nst;
 463         }
 464
 465         if (ir->epc != epcNO)
 466         {
 467             if (ir->nstpcouple < 0)
 468             {
 469                 ir->nstpcouple = ir_optimal_nstpcouple(ir);
 470             }
 471         }
 472         if (IR_TWINRANGE(*ir))
 473         {
 474             check_nst("nstlist", ir->nstlist,
 475                       "nstcalcenergy", &ir->nstcalcenergy, wi);
 476             if (ir->epc != epcNO)
 477             {
 478                 check_nst("nstlist", ir->nstlist,
 479                           "nstpcouple", &ir->nstpcouple, wi);
 480             }
 481         }
 482
 483         if (ir->nstcalcenergy > 0)
 484         {
 485             if (ir->efep != efepNO)
 486             {
 487                 /* nstdhdl should be a multiple of nstcalcenergy */
 488                 check_nst("nstcalcenergy", ir->nstcalcenergy,
 489                           "nstdhdl", &ir->fepvals->nstdhdl, wi);
 490                 /* nstexpanded should be a multiple of nstcalcenergy */
 491                 check_nst("nstcalcenergy", ir->nstcalcenergy,
 492                           "nstexpanded", &ir->expandedvals->nstexpanded, wi);
 493             }
 494             /* for storing exact averages nstenergy should be
 495              * a multiple of nstcalcenergy
 496              */
 497             check_nst("nstcalcenergy", ir->nstcalcenergy,
 498                       "nstenergy", &ir->nstenergy, wi);
 499         }
 500     }
 501
 502     /* LD STUFF */
 503     if ((EI_SD(ir->eI) || ir->eI == eiBD) &&
 504         ir->bContinuation && ir->ld_seed != -1)
 505     {
 506         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)");
 507     }
 508
 509     /* TPI STUFF */
 510     if (EI_TPI(ir->eI))
 511     {
 512         sprintf(err_buf, "TPI only works with pbc = %s", epbc_names[epbcXYZ]);
 513         CHECK(ir->ePBC != epbcXYZ);
 514         sprintf(err_buf, "TPI only works with ns = %s", ens_names[ensGRID]);
 515         CHECK(ir->ns_type != ensGRID);
 516         sprintf(err_buf, "with TPI nstlist should be larger than zero");
 517         CHECK(ir->nstlist <= 0);
 518         sprintf(err_buf, "TPI does not work with full electrostatics other than PME");
 519         CHECK(EEL_FULL(ir->coulombtype) && !EEL_PME(ir->coulombtype));
 520     }
 521
 522     /* SHAKE / LINCS */
 523     if ( (opts->nshake > 0) && (opts->bMorse) )
 524     {
 525         sprintf(warn_buf,
 526                 "Using morse bond-potentials while constraining bonds is useless");
 527         warning(wi, warn_buf);
 528     }
 529
 530     if ((EI_SD(ir->eI) || ir->eI == eiBD) &&
 531         ir->bContinuation && ir->ld_seed != -1)
 532     {
 533         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)");
 534     }
 535     /* verify simulated tempering options */
 536
 537     if (ir->bSimTemp)
 538     {
 539         gmx_bool bAllTempZero = TRUE;
 540         for (i = 0; i < fep->n_lambda; i++)
 541         {
 542             sprintf(err_buf, "Entry %d for %s must be between 0 and 1, instead is %g", i, efpt_names[efptTEMPERATURE], fep->all_lambda[efptTEMPERATURE][i]);
 543             CHECK((fep->all_lambda[efptTEMPERATURE][i] < 0) || (fep->all_lambda[efptTEMPERATURE][i] > 1));
 544             if (fep->all_lambda[efptTEMPERATURE][i] > 0)
 545             {
 546                 bAllTempZero = FALSE;
 547             }
 548         }
 549         sprintf(err_buf, "if simulated tempering is on, temperature-lambdas may not be all zero");
 550         CHECK(bAllTempZero == TRUE);
 551
 552         sprintf(err_buf, "Simulated tempering is currently only compatible with md-vv");
 553         CHECK(ir->eI != eiVV);
 554
 555         /* check compatability of the temperature coupling with simulated tempering */
 556
 557         if (ir->etc == etcNOSEHOOVER)
 558         {
 559             sprintf(warn_buf, "Nose-Hoover based temperature control such as [%s] my not be entirelyconsistent with simulated tempering", etcoupl_names[ir->etc]);
 560             warning_note(wi, warn_buf);
 561         }
 562
 563         /* check that the temperatures make sense */
 564
 565         sprintf(err_buf, "Higher simulated tempering temperature (%g) must be >= than the simulated tempering lower temperature (%g)", ir->simtempvals->simtemp_high, ir->simtempvals->simtemp_low);
 566         CHECK(ir->simtempvals->simtemp_high <= ir->simtempvals->simtemp_low);
 567
 568         sprintf(err_buf, "Higher simulated tempering temperature (%g) must be >= zero", ir->simtempvals->simtemp_high);
 569         CHECK(ir->simtempvals->simtemp_high <= 0);
 570
 571         sprintf(err_buf, "Lower simulated tempering temperature (%g) must be >= zero", ir->simtempvals->simtemp_low);
 572         CHECK(ir->simtempvals->simtemp_low <= 0);
 573     }
 574
 575     /* verify free energy options */
 576
 577     if (ir->efep != efepNO)
 578     {
 579         fep = ir->fepvals;
 580         sprintf(err_buf, "The soft-core power is %d and can only be 1 or 2",
 581                 fep->sc_power);
 582         CHECK(fep->sc_alpha != 0 && fep->sc_power != 1 && fep->sc_power != 2);
 583
 584         sprintf(err_buf, "The soft-core sc-r-power is %d and can only be 6 or 48",
 585                 (int)fep->sc_r_power);
 586         CHECK(fep->sc_alpha != 0 && fep->sc_r_power != 6.0 && fep->sc_r_power != 48.0);
 587
 588         sprintf(err_buf, "Can't use postive delta-lambda (%g) if initial state/lambda does not start at zero", fep->delta_lambda);
 589         CHECK(fep->delta_lambda > 0 && ((fep->init_fep_state > 0) ||  (fep->init_lambda > 0)));
 590
 591         sprintf(err_buf, "Can't use postive delta-lambda (%g) with expanded ensemble simulations", fep->delta_lambda);
 592         CHECK(fep->delta_lambda > 0 && (ir->efep == efepEXPANDED));
 593
 594         sprintf(err_buf, "Free-energy not implemented for Ewald");
 595         CHECK(ir->coulombtype == eelEWALD);
 596
 597         /* check validty of lambda inputs */
 598         if (fep->n_lambda == 0)
 599         {
 600             /* Clear output in case of no states:*/
 601             sprintf(err_buf, "init-lambda-state set to %d: no lambda states are defined.", fep->init_fep_state);
 602             CHECK((fep->init_fep_state >= 0) && (fep->n_lambda == 0));
 603         }
 604         else
 605         {
 606             sprintf(err_buf, "initial thermodynamic state %d does not exist, only goes to %d", fep->init_fep_state, fep->n_lambda-1);
 607             CHECK((fep->init_fep_state >= fep->n_lambda));
 608         }
 609
 610         sprintf(err_buf, "Lambda state must be set, either with init-lambda-state or with init-lambda");
 611         CHECK((fep->init_fep_state < 0) && (fep->init_lambda < 0));
 612
 613         sprintf(err_buf, "init-lambda=%g while init-lambda-state=%d. Lambda state must be set either with init-lambda-state or with init-lambda, but not both",
 614                 fep->init_lambda, fep->init_fep_state);
 615         CHECK((fep->init_fep_state >= 0) && (fep->init_lambda >= 0));
 616
 617
 618
 619         if ((fep->init_lambda >= 0) && (fep->delta_lambda == 0))
 620         {
 621             int n_lambda_terms;
 622             n_lambda_terms = 0;
 623             for (i = 0; i < efptNR; i++)
 624             {
 625                 if (fep->separate_dvdl[i])
 626                 {
 627                     n_lambda_terms++;
 628                 }
 629             }
 630             if (n_lambda_terms > 1)
 631             {
 632                 sprintf(warn_buf, "If lambda vector states (fep-lambdas, coul-lambdas etc.) are set, don't use init-lambda to set lambda state (except for slow growth). Use init-lambda-state instead.");
 633                 warning(wi, warn_buf);
 634             }
 635
 636             if (n_lambda_terms < 2 && fep->n_lambda > 0)
 637             {
 638                 warning_note(wi,
 639                              "init-lambda is deprecated for setting lambda state (except for slow growth). Use init-lambda-state instead.");
 640             }
 641         }
 642
 643         for (j = 0; j < efptNR; j++)
 644         {
 645             for (i = 0; i < fep->n_lambda; i++)
 646             {
 647                 sprintf(err_buf, "Entry %d for %s must be between 0 and 1, instead is %g", i, efpt_names[j], fep->all_lambda[j][i]);
 648                 CHECK((fep->all_lambda[j][i] < 0) || (fep->all_lambda[j][i] > 1));
 649             }
 650         }
 651
 652         if ((fep->sc_alpha > 0) && (!fep->bScCoul))
 653         {
 654             for (i = 0; i < fep->n_lambda; i++)
 655             {
 656                 sprintf(err_buf, "For state %d, vdw-lambdas (%f) is changing with vdw softcore, while coul-lambdas (%f) is nonzero without coulomb softcore: this will lead to crashes, and is not supported.", i, fep->all_lambda[efptVDW][i],
 657                         fep->all_lambda[efptCOUL][i]);
 658                 CHECK((fep->sc_alpha > 0) &&
 659                       (((fep->all_lambda[efptCOUL][i] > 0.0) &&
 660                         (fep->all_lambda[efptCOUL][i] < 1.0)) &&
 661                        ((fep->all_lambda[efptVDW][i] > 0.0) &&
 662                         (fep->all_lambda[efptVDW][i] < 1.0))));
 663             }
 664         }
 665
 666         if ((fep->bScCoul) && (EEL_PME(ir->coulombtype)))
 667         {
 668             real sigma, lambda, r_sc;
 669
 670             sigma  = 0.34;
 671             /* Maximum estimate for A and B charges equal with lambda power 1 */
 672             lambda = 0.5;
 673             r_sc   = pow(lambda*fep->sc_alpha*pow(sigma/ir->rcoulomb, fep->sc_r_power) + 1.0, 1.0/fep->sc_r_power);
 674             sprintf(warn_buf, "With PME there is a minor soft core effect present at the cut-off, proportional to (LJsigma/rcoulomb)^%g. This could have a minor effect on energy conservation, but usually other effects dominate. With a common sigma value of %g nm the fraction of the particle-particle potential at the cut-off at lambda=%g is around %.1e, while ewald-rtol is %.1e.",
 675                     fep->sc_r_power,
 676                     sigma, lambda, r_sc - 1.0, ir->ewald_rtol);
 677             warning_note(wi, warn_buf);
 678         }
 679
 680         /*  Free Energy Checks -- In an ideal world, slow growth and FEP would
 681             be treated differently, but that's the next step */
 682
 683         for (i = 0; i < efptNR; i++)
 684         {
 685             for (j = 0; j < fep->n_lambda; j++)
 686             {
 687                 sprintf(err_buf, "%s[%d] must be between 0 and 1", efpt_names[i], j);
 688                 CHECK((fep->all_lambda[i][j] < 0) || (fep->all_lambda[i][j] > 1));
 689             }
 690         }
 691     }
 692
 693     if ((ir->bSimTemp) || (ir->efep == efepEXPANDED))
 694     {
 695         fep    = ir->fepvals;
 696         expand = ir->expandedvals;
 697
 698         /* checking equilibration of weights inputs for validity */
 699
 700         sprintf(err_buf, "weight-equil-number-all-lambda (%d) is ignored if lmc-weights-equil is not equal to %s",
 701                 expand->equil_n_at_lam, elmceq_names[elmceqNUMATLAM]);
 702         CHECK((expand->equil_n_at_lam > 0) && (expand->elmceq != elmceqNUMATLAM));
 703
 704         sprintf(err_buf, "weight-equil-number-samples (%d) is ignored if lmc-weights-equil is not equal to %s",
 705                 expand->equil_samples, elmceq_names[elmceqSAMPLES]);
 706         CHECK((expand->equil_samples > 0) && (expand->elmceq != elmceqSAMPLES));
 707
 708         sprintf(err_buf, "weight-equil-number-steps (%d) is ignored if lmc-weights-equil is not equal to %s",
 709                 expand->equil_steps, elmceq_names[elmceqSTEPS]);
 710         CHECK((expand->equil_steps > 0) && (expand->elmceq != elmceqSTEPS));
 711
 712         sprintf(err_buf, "weight-equil-wl-delta (%d) is ignored if lmc-weights-equil is not equal to %s",
 713                 expand->equil_samples, elmceq_names[elmceqWLDELTA]);
 714         CHECK((expand->equil_wl_delta > 0) && (expand->elmceq != elmceqWLDELTA));
 715
 716         sprintf(err_buf, "weight-equil-count-ratio (%f) is ignored if lmc-weights-equil is not equal to %s",
 717                 expand->equil_ratio, elmceq_names[elmceqRATIO]);
 718         CHECK((expand->equil_ratio > 0) && (expand->elmceq != elmceqRATIO));
 719
 720         sprintf(err_buf, "weight-equil-number-all-lambda (%d) must be a positive integer if lmc-weights-equil=%s",
 721                 expand->equil_n_at_lam, elmceq_names[elmceqNUMATLAM]);
 722         CHECK((expand->equil_n_at_lam <= 0) && (expand->elmceq == elmceqNUMATLAM));
 723
 724         sprintf(err_buf, "weight-equil-number-samples (%d) must be a positive integer if lmc-weights-equil=%s",
 725                 expand->equil_samples, elmceq_names[elmceqSAMPLES]);
 726         CHECK((expand->equil_samples <= 0) && (expand->elmceq == elmceqSAMPLES));
 727
 728         sprintf(err_buf, "weight-equil-number-steps (%d) must be a positive integer if lmc-weights-equil=%s",
 729                 expand->equil_steps, elmceq_names[elmceqSTEPS]);
 730         CHECK((expand->equil_steps <= 0) && (expand->elmceq == elmceqSTEPS));
 731
 732         sprintf(err_buf, "weight-equil-wl-delta (%f) must be > 0 if lmc-weights-equil=%s",
 733                 expand->equil_wl_delta, elmceq_names[elmceqWLDELTA]);
 734         CHECK((expand->equil_wl_delta <= 0) && (expand->elmceq == elmceqWLDELTA));
 735
 736         sprintf(err_buf, "weight-equil-count-ratio (%f) must be > 0 if lmc-weights-equil=%s",
 737                 expand->equil_ratio, elmceq_names[elmceqRATIO]);
 738         CHECK((expand->equil_ratio <= 0) && (expand->elmceq == elmceqRATIO));
 739
 740         sprintf(err_buf, "lmc-weights-equil=%s only possible when lmc-stats = %s or lmc-stats %s",
 741                 elmceq_names[elmceqWLDELTA], elamstats_names[elamstatsWL], elamstats_names[elamstatsWWL]);
 742         CHECK((expand->elmceq == elmceqWLDELTA) && (!EWL(expand->elamstats)));
 743
 744         sprintf(err_buf, "lmc-repeats (%d) must be greater than 0", expand->lmc_repeats);
 745         CHECK((expand->lmc_repeats <= 0));
 746         sprintf(err_buf, "minimum-var-min (%d) must be greater than 0", expand->minvarmin);
 747         CHECK((expand->minvarmin <= 0));
 748         sprintf(err_buf, "weight-c-range (%d) must be greater or equal to 0", expand->c_range);
 749         CHECK((expand->c_range < 0));
 750         sprintf(err_buf, "init-lambda-state (%d) must be zero if lmc-forced-nstart (%d)> 0 and lmc-move != 'no'",
 751                 fep->init_fep_state, expand->lmc_forced_nstart);
 752         CHECK((fep->init_fep_state != 0) && (expand->lmc_forced_nstart > 0) && (expand->elmcmove != elmcmoveNO));
 753         sprintf(err_buf, "lmc-forced-nstart (%d) must not be negative", expand->lmc_forced_nstart);
 754         CHECK((expand->lmc_forced_nstart < 0));
 755         sprintf(err_buf, "init-lambda-state (%d) must be in the interval [0,number of lambdas)", fep->init_fep_state);
 756         CHECK((fep->init_fep_state < 0) || (fep->init_fep_state >= fep->n_lambda));
 757
 758         sprintf(err_buf, "init-wl-delta (%f) must be greater than or equal to 0", expand->init_wl_delta);
 759         CHECK((expand->init_wl_delta < 0));
 760         sprintf(err_buf, "wl-ratio (%f) must be between 0 and 1", expand->wl_ratio);
 761         CHECK((expand->wl_ratio <= 0) || (expand->wl_ratio >= 1));
 762         sprintf(err_buf, "wl-scale (%f) must be between 0 and 1", expand->wl_scale);
 763         CHECK((expand->wl_scale <= 0) || (expand->wl_scale >= 1));
 764
 765         /* if there is no temperature control, we need to specify an MC temperature */
 766         sprintf(err_buf, "If there is no temperature control, and lmc-mcmove!= 'no',mc_temperature must be set to a positive number");
 767         if (expand->nstTij > 0)
 768         {
 769             sprintf(err_buf, "nst-transition-matrix (%d) must be an integer multiple of nstlog (%d)",
 770                     expand->nstTij, ir->nstlog);
 771             CHECK((mod(expand->nstTij, ir->nstlog) != 0));
 772         }
 773     }
 774
 775     /* PBC/WALLS */
 776     sprintf(err_buf, "walls only work with pbc=%s", epbc_names[epbcXY]);
 777     CHECK(ir->nwall && ir->ePBC != epbcXY);
 778
 779     /* VACUUM STUFF */
 780     if (ir->ePBC != epbcXYZ && ir->nwall != 2)
 781     {
 782         if (ir->ePBC == epbcNONE)
 783         {
 784             if (ir->epc != epcNO)
 785             {
 786                 warning(wi, "Turning off pressure coupling for vacuum system");
 787                 ir->epc = epcNO;
 788             }
 789         }
 790         else
 791         {
 792             sprintf(err_buf, "Can not have pressure coupling with pbc=%s",
 793                     epbc_names[ir->ePBC]);
 794             CHECK(ir->epc != epcNO);
 795         }
 796         sprintf(err_buf, "Can not have Ewald with pbc=%s", epbc_names[ir->ePBC]);
 797         CHECK(EEL_FULL(ir->coulombtype));
 798
 799         sprintf(err_buf, "Can not have dispersion correction with pbc=%s",
 800                 epbc_names[ir->ePBC]);
 801         CHECK(ir->eDispCorr != edispcNO);
 802     }
 803
 804     if (ir->rlist == 0.0)
 805     {
 806         sprintf(err_buf, "can only have neighborlist cut-off zero (=infinite)\n"
 807                 "with coulombtype = %s or coulombtype = %s\n"
 808                 "without periodic boundary conditions (pbc = %s) and\n"
 809                 "rcoulomb and rvdw set to zero",
 810                 eel_names[eelCUT], eel_names[eelUSER], epbc_names[epbcNONE]);
 811         CHECK(((ir->coulombtype != eelCUT) && (ir->coulombtype != eelUSER)) ||
 812               (ir->ePBC     != epbcNONE) ||
 813               (ir->rcoulomb != 0.0)      || (ir->rvdw != 0.0));
 814
 815         if (ir->nstlist < 0)
 816         {
 817             warning_error(wi, "Can not have heuristic neighborlist updates without cut-off");
 818         }
 819         if (ir->nstlist > 0)
 820         {
 821             warning_note(wi, "Simulating without cut-offs is usually (slightly) faster with nstlist=0, nstype=simple and particle decomposition");
 822         }
 823     }
 824
 825     /* COMM STUFF */
 826     if (ir->nstcomm == 0)
 827     {
 828         ir->comm_mode = ecmNO;
 829     }
 830     if (ir->comm_mode != ecmNO)
 831     {
 832         if (ir->nstcomm < 0)
 833         {
 834             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");
 835             ir->nstcomm = abs(ir->nstcomm);
 836         }
 837
 838         if (ir->nstcalcenergy > 0 && ir->nstcomm < ir->nstcalcenergy)
 839         {
 840             warning_note(wi, "nstcomm < nstcalcenergy defeats the purpose of nstcalcenergy, setting nstcomm to nstcalcenergy");
 841             ir->nstcomm = ir->nstcalcenergy;
 842         }
 843
 844         if (ir->comm_mode == ecmANGULAR)
 845         {
 846             sprintf(err_buf, "Can not remove the rotation around the center of mass with periodic molecules");
 847             CHECK(ir->bPeriodicMols);
 848             if (ir->ePBC != epbcNONE)
 849             {
 850                 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).");
 851             }
 852         }
 853     }
 854
 855     if (EI_STATE_VELOCITY(ir->eI) && ir->ePBC == epbcNONE && ir->comm_mode != ecmANGULAR)
 856     {
 857         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.");
 858     }
 859
 860     sprintf(err_buf, "Twin-range neighbour searching (NS) with simple NS"
 861             " algorithm not implemented");
 862     CHECK(((ir->rcoulomb > ir->rlist) || (ir->rvdw > ir->rlist))
 863           && (ir->ns_type == ensSIMPLE));
 864
 865     /* TEMPERATURE COUPLING */
 866     if (ir->etc == etcYES)
 867     {
 868         ir->etc = etcBERENDSEN;
 869         warning_note(wi, "Old option for temperature coupling given: "
 870                      "changing \"yes\" to \"Berendsen\"\n");
 871     }
 872
 873     if ((ir->etc == etcNOSEHOOVER) || (ir->epc == epcMTTK))
 874     {
 875         if (ir->opts.nhchainlength < 1)
 876         {
 877             sprintf(warn_buf, "number of Nose-Hoover chains (currently %d) cannot be less than 1,reset to 1\n", ir->opts.nhchainlength);
 878             ir->opts.nhchainlength = 1;
 879             warning(wi, warn_buf);
 880         }
 881
 882         if (ir->etc == etcNOSEHOOVER && !EI_VV(ir->eI) && ir->opts.nhchainlength > 1)
 883         {
 884             warning_note(wi, "leapfrog does not yet support Nose-Hoover chains, nhchainlength reset to 1");
 885             ir->opts.nhchainlength = 1;
 886         }
 887     }
 888     else
 889     {
 890         ir->opts.nhchainlength = 0;
 891     }
 892
 893     if (ir->eI == eiVVAK)
 894     {
 895         sprintf(err_buf, "%s implemented primarily for validation, and requires nsttcouple = 1 and nstpcouple = 1.",
 896                 ei_names[eiVVAK]);
 897         CHECK((ir->nsttcouple != 1) || (ir->nstpcouple != 1));
 898     }
 899
 900     if (ETC_ANDERSEN(ir->etc))
 901     {
 902         sprintf(err_buf, "%s temperature control not supported for integrator %s.", etcoupl_names[ir->etc], ei_names[ir->eI]);
 903         CHECK(!(EI_VV(ir->eI)));
 904
 905         for (i = 0; i < ir->opts.ngtc; i++)
 906         {
 907             sprintf(err_buf, "all tau_t must currently be equal using Andersen temperature control, violated for group %d", i);
 908             CHECK(ir->opts.tau_t[0] != ir->opts.tau_t[i]);
 909             sprintf(err_buf, "all tau_t must be postive using Andersen temperature control, tau_t[%d]=%10.6f",
 910                     i, ir->opts.tau_t[i]);
 911             CHECK(ir->opts.tau_t[i] < 0);
 912         }
 913         if (ir->nstcomm > 0 && (ir->etc == etcANDERSEN))
 914         {
 915             sprintf(warn_buf, "Center of mass removal not necessary for %s.  All velocities of coupled groups are rerandomized periodically, so flying ice cube errors will not occur.", etcoupl_names[ir->etc]);
 916             warning_note(wi, warn_buf);
 917         }
 918
 919         sprintf(err_buf, "nstcomm must be 1, not %d for %s, as velocities of atoms in coupled groups are randomized every time step", ir->nstcomm, etcoupl_names[ir->etc]);
 920         CHECK(ir->nstcomm > 1 && (ir->etc == etcANDERSEN));
 921
 922         for (i = 0; i < ir->opts.ngtc; i++)
 923         {
 924             int nsteps = (int)(ir->opts.tau_t[i]/ir->delta_t);
 925             sprintf(err_buf, "tau_t/delta_t for group %d for temperature control method %s must be a multiple of nstcomm (%d), as velocities of atoms in coupled groups are randomized every time step. The input tau_t (%8.3f) leads to %d steps per randomization", i, etcoupl_names[ir->etc], ir->nstcomm, ir->opts.tau_t[i], nsteps);
 926             CHECK((nsteps % ir->nstcomm) && (ir->etc == etcANDERSENMASSIVE));
 927         }
 928     }
 929     if (ir->etc == etcBERENDSEN)
 930     {
 931         sprintf(warn_buf, "The %s thermostat does not generate the correct kinetic energy distribution. You might want to consider using the %s thermostat.",
 932                 ETCOUPLTYPE(ir->etc), ETCOUPLTYPE(etcVRESCALE));
 933         warning_note(wi, warn_buf);
 934     }
 935
 936     if ((ir->etc == etcNOSEHOOVER || ETC_ANDERSEN(ir->etc))
 937         && ir->epc == epcBERENDSEN)
 938     {
 939         sprintf(warn_buf, "Using Berendsen pressure coupling invalidates the "
 940                 "true ensemble for the thermostat");
 941         warning(wi, warn_buf);
 942     }
 943
 944     /* PRESSURE COUPLING */
 945     if (ir->epc == epcISOTROPIC)
 946     {
 947         ir->epc = epcBERENDSEN;
 948         warning_note(wi, "Old option for pressure coupling given: "
 949                      "changing \"Isotropic\" to \"Berendsen\"\n");
 950     }
 951
 952     if (ir->epc != epcNO)
 953     {
 954         dt_pcoupl = ir->nstpcouple*ir->delta_t;
 955
 956         sprintf(err_buf, "tau-p must be > 0 instead of %g\n", ir->tau_p);
 957         CHECK(ir->tau_p <= 0);
 958
 959         if (ir->tau_p/dt_pcoupl < pcouple_min_integration_steps(ir->epc))
 960         {
 961             sprintf(warn_buf, "For proper integration of the %s barostat, tau-p (%g) should be at least %d times larger than nstpcouple*dt (%g)",
 962                     EPCOUPLTYPE(ir->epc), ir->tau_p, pcouple_min_integration_steps(ir->epc), dt_pcoupl);
 963             warning(wi, warn_buf);
 964         }
 965
 966         sprintf(err_buf, "compressibility must be > 0 when using pressure"
 967                 " coupling %s\n", EPCOUPLTYPE(ir->epc));
 968         CHECK(ir->compress[XX][XX] < 0 || ir->compress[YY][YY] < 0 ||
 969               ir->compress[ZZ][ZZ] < 0 ||
 970               (trace(ir->compress) == 0 && ir->compress[YY][XX] <= 0 &&
 971                ir->compress[ZZ][XX] <= 0 && ir->compress[ZZ][YY] <= 0));
 972
 973         if (epcPARRINELLORAHMAN == ir->epc && opts->bGenVel)
 974         {
 975             sprintf(warn_buf,
 976                     "You are generating velocities so I am assuming you "
 977                     "are equilibrating a system. You are using "
 978                     "%s pressure coupling, but this can be "
 979                     "unstable for equilibration. If your system crashes, try "
 980                     "equilibrating first with Berendsen pressure coupling. If "
 981                     "you are not equilibrating the system, you can probably "
 982                     "ignore this warning.",
 983                     epcoupl_names[ir->epc]);
 984             warning(wi, warn_buf);
 985         }
 986     }
 987
 988     if (EI_VV(ir->eI))
 989     {
 990         if (ir->epc > epcNO)
 991         {
 992             if ((ir->epc != epcBERENDSEN) && (ir->epc != epcMTTK))
 993             {
 994                 warning_error(wi, "for md-vv and md-vv-avek, can only use Berendsen and Martyna-Tuckerman-Tobias-Klein (MTTK) equations for pressure control; MTTK is equivalent to Parrinello-Rahman.");
 995             }
 996         }
 997     }
 998
 999     /* ELECTROSTATICS */
1000     /* More checks are in triple check (grompp.c) */
1001
1002     if (ir->coulombtype == eelSWITCH)
1003     {
1004         sprintf(warn_buf, "coulombtype = %s is only for testing purposes and can lead to serious "
1005                 "artifacts, advice: use coulombtype = %s",
1006                 eel_names[ir->coulombtype],
1007                 eel_names[eelRF_ZERO]);
1008         warning(wi, warn_buf);
1009     }
1010
1011     if (ir->epsilon_r != 1 && ir->implicit_solvent == eisGBSA)
1012     {
1013         sprintf(warn_buf, "epsilon-r = %g with GB implicit solvent, will use this value for inner dielectric", ir->epsilon_r);
1014         warning_note(wi, warn_buf);
1015     }
1016
1017     if (EEL_RF(ir->coulombtype) && ir->epsilon_rf == 1 && ir->epsilon_r != 1)
1018     {
1019         sprintf(warn_buf, "epsilon-r = %g and epsilon-rf = 1 with reaction field, proceeding assuming old format and exchanging epsilon-r and epsilon-rf", ir->epsilon_r);
1020         warning(wi, warn_buf);
1021         ir->epsilon_rf = ir->epsilon_r;
1022         ir->epsilon_r  = 1.0;
1023     }
1024
1025     if (getenv("GALACTIC_DYNAMICS") == NULL)
1026     {
1027         sprintf(err_buf, "epsilon-r must be >= 0 instead of %g\n", ir->epsilon_r);
1028         CHECK(ir->epsilon_r < 0);
1029     }
1030
1031     if (EEL_RF(ir->coulombtype))
1032     {
1033         /* reaction field (at the cut-off) */
1034
1035         if (ir->coulombtype == eelRF_ZERO)
1036         {
1037             sprintf(warn_buf, "With coulombtype = %s, epsilon-rf must be 0, assuming you meant epsilon_rf=0",
1038                     eel_names[ir->coulombtype]);
1039             CHECK(ir->epsilon_rf != 0);
1040             ir->epsilon_rf = 0.0;
1041         }
1042
1043         sprintf(err_buf, "epsilon-rf must be >= epsilon-r");
1044         CHECK((ir->epsilon_rf < ir->epsilon_r && ir->epsilon_rf != 0) ||
1045               (ir->epsilon_r == 0));
1046         if (ir->epsilon_rf == ir->epsilon_r)
1047         {
1048             sprintf(warn_buf, "Using epsilon-rf = epsilon-r with %s does not make sense",
1049                     eel_names[ir->coulombtype]);
1050             warning(wi, warn_buf);
1051         }
1052     }
1053     /* Allow rlist>rcoulomb for tabulated long range stuff. This just
1054      * means the interaction is zero outside rcoulomb, but it helps to
1055      * provide accurate energy conservation.
1056      */
1057     if (EEL_MIGHT_BE_ZERO_AT_CUTOFF(ir->coulombtype))
1058     {
1059         if (EEL_SWITCHED(ir->coulombtype))
1060         {
1061             sprintf(err_buf,
1062                     "With coulombtype = %s rcoulomb_switch must be < rcoulomb. Or, better: Use the potential modifier options!",
1063                     eel_names[ir->coulombtype]);
1064             CHECK(ir->rcoulomb_switch >= ir->rcoulomb);
1065         }
1066     }
1067     else if (ir->coulombtype == eelCUT || EEL_RF(ir->coulombtype))
1068     {
1069         if (ir->cutoff_scheme == ecutsGROUP && ir->coulomb_modifier == eintmodNONE)
1070         {
1071             sprintf(err_buf, "With coulombtype = %s, rcoulomb should be >= rlist unless you use a potential modifier",
1072                     eel_names[ir->coulombtype]);
1073             CHECK(ir->rlist > ir->rcoulomb);
1074         }
1075     }
1076
1077     if (ir->coulombtype == eelSWITCH || ir->coulombtype == eelSHIFT ||
1078         ir->vdwtype == evdwSWITCH || ir->vdwtype == evdwSHIFT)
1079     {
1080         sprintf(warn_buf,
1081                 "The switch/shift interaction settings are just for compatibility; you will get better "
1082                 "performance from applying potential modifiers to your interactions!\n");
1083         warning_note(wi, warn_buf);
1084     }
1085
1086     if (ir->coulombtype == eelPMESWITCH)
1087     {
1088         if (ir->rcoulomb_switch/ir->rcoulomb < 0.9499)
1089         {
1090             sprintf(warn_buf, "The switching range for %s should be 5%% or less, energy conservation will be good anyhow, since ewald_rtol = %g",
1091                     eel_names[ir->coulombtype],
1092                     ir->ewald_rtol);
1093             warning(wi, warn_buf);
1094         }
1095     }
1096
1097     if (EEL_FULL(ir->coulombtype))
1098     {
1099         if (ir->coulombtype == eelPMESWITCH || ir->coulombtype == eelPMEUSER ||
1100             ir->coulombtype == eelPMEUSERSWITCH)
1101         {
1102             sprintf(err_buf, "With coulombtype = %s, rcoulomb must be <= rlist",
1103                     eel_names[ir->coulombtype]);
1104             CHECK(ir->rcoulomb > ir->rlist);
1105         }
1106         else if (ir->cutoff_scheme == ecutsGROUP && ir->coulomb_modifier == eintmodNONE)
1107         {
1108             if (ir->coulombtype == eelPME || ir->coulombtype == eelP3M_AD)
1109             {
1110                 sprintf(err_buf,
1111                         "With coulombtype = %s (without modifier), rcoulomb must be equal to rlist,\n"
1112                         "or rlistlong if nstcalclr=1. For optimal energy conservation,consider using\n"
1113                         "a potential modifier.", eel_names[ir->coulombtype]);
1114                 if (ir->nstcalclr == 1)
1115                 {
1116                     CHECK(ir->rcoulomb != ir->rlist && ir->rcoulomb != ir->rlistlong);
1117                 }
1118                 else
1119                 {
1120                     CHECK(ir->rcoulomb != ir->rlist);
1121                 }
1122             }
1123         }
1124     }
1125
1126     if (EEL_PME(ir->coulombtype))
1127     {
1128         if (ir->pme_order < 3)
1129         {
1130             warning_error(wi, "pme-order can not be smaller than 3");
1131         }
1132     }
1133
1134     if (ir->nwall == 2 && EEL_FULL(ir->coulombtype))
1135     {
1136         if (ir->ewald_geometry == eewg3D)
1137         {
1138             sprintf(warn_buf, "With pbc=%s you should use ewald-geometry=%s",
1139                     epbc_names[ir->ePBC], eewg_names[eewg3DC]);
1140             warning(wi, warn_buf);
1141         }
1142         /* This check avoids extra pbc coding for exclusion corrections */
1143         sprintf(err_buf, "wall-ewald-zfac should be >= 2");
1144         CHECK(ir->wall_ewald_zfac < 2);
1145     }
1146
1147     if (EVDW_SWITCHED(ir->vdwtype))
1148     {
1149         sprintf(err_buf, "With vdwtype = %s rvdw-switch must be < rvdw. Or, better - use a potential modifier.",
1150                 evdw_names[ir->vdwtype]);
1151         CHECK(ir->rvdw_switch >= ir->rvdw);
1152     }
1153     else if (ir->vdwtype == evdwCUT)
1154     {
1155         if (ir->cutoff_scheme == ecutsGROUP && ir->vdw_modifier == eintmodNONE)
1156         {
1157             sprintf(err_buf, "With vdwtype = %s, rvdw must be >= rlist unless you use a potential modifier", evdw_names[ir->vdwtype]);
1158             CHECK(ir->rlist > ir->rvdw);
1159         }
1160     }
1161     if (ir->cutoff_scheme == ecutsGROUP)
1162     {
1163         if (((ir->coulomb_modifier != eintmodNONE && ir->rcoulomb == ir->rlist) ||
1164              (ir->vdw_modifier != eintmodNONE && ir->rvdw == ir->rlist)) &&
1165             ir->nstlist != 1)
1166         {
1167             warning_note(wi, "With exact cut-offs, rlist should be "
1168                          "larger than rcoulomb and rvdw, so that there "
1169                          "is a buffer region for particle motion "
1170                          "between neighborsearch steps");
1171         }
1172
1173         if (EEL_IS_ZERO_AT_CUTOFF(ir->coulombtype)
1174             && (ir->rlistlong <= ir->rcoulomb))
1175         {
1176             sprintf(warn_buf, "For energy conservation with switch/shift potentials, %s should be 0.1 to 0.3 nm larger than rcoulomb.",
1177                     IR_TWINRANGE(*ir) ? "rlistlong" : "rlist");
1178             warning_note(wi, warn_buf);
1179         }
1180         if (EVDW_SWITCHED(ir->vdwtype) && (ir->rlistlong <= ir->rvdw))
1181         {
1182             sprintf(warn_buf, "For energy conservation with switch/shift potentials, %s should be 0.1 to 0.3 nm larger than rvdw.",
1183                     IR_TWINRANGE(*ir) ? "rlistlong" : "rlist");
1184             warning_note(wi, warn_buf);
1185         }
1186     }
1187
1188     if (ir->vdwtype == evdwUSER && ir->eDispCorr != edispcNO)
1189     {
1190         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.");
1191     }
1192
1193     if (ir->nstlist == -1)
1194     {
1195         sprintf(err_buf, "With nstlist=-1 rvdw and rcoulomb should be smaller than rlist to account for diffusion and possibly charge-group radii");
1196         CHECK(ir->rvdw >= ir->rlist || ir->rcoulomb >= ir->rlist);
1197     }
1198     sprintf(err_buf, "nstlist can not be smaller than -1");
1199     CHECK(ir->nstlist < -1);
1200
1201     if (ir->eI == eiLBFGS && (ir->coulombtype == eelCUT || ir->vdwtype == evdwCUT)
1202         && ir->rvdw != 0)
1203     {
1204         warning(wi, "For efficient BFGS minimization, use switch/shift/pme instead of cut-off.");
1205     }
1206
1207     if (ir->eI == eiLBFGS && ir->nbfgscorr <= 0)
1208     {
1209         warning(wi, "Using L-BFGS with nbfgscorr<=0 just gets you steepest descent.");
1210     }
1211
1212     /* ENERGY CONSERVATION */
1213     if (ir_NVE(ir) && ir->cutoff_scheme == ecutsGROUP)
1214     {
1215         if (!EVDW_MIGHT_BE_ZERO_AT_CUTOFF(ir->vdwtype) && ir->rvdw > 0 && ir->vdw_modifier == eintmodNONE)
1216         {
1217             sprintf(warn_buf, "You are using a cut-off for VdW interactions with NVE, for good energy conservation use vdwtype = %s (possibly with DispCorr)",
1218                     evdw_names[evdwSHIFT]);
1219             warning_note(wi, warn_buf);
1220         }
1221         if (!EEL_MIGHT_BE_ZERO_AT_CUTOFF(ir->coulombtype) && ir->rcoulomb > 0 && ir->coulomb_modifier == eintmodNONE)
1222         {
1223             sprintf(warn_buf, "You are using a cut-off for electrostatics with NVE, for good energy conservation use coulombtype = %s or %s",
1224                     eel_names[eelPMESWITCH], eel_names[eelRF_ZERO]);
1225             warning_note(wi, warn_buf);
1226         }
1227     }
1228
1229     /* IMPLICIT SOLVENT */
1230     if (ir->coulombtype == eelGB_NOTUSED)
1231     {
1232         ir->coulombtype      = eelCUT;
1233         ir->implicit_solvent = eisGBSA;
1234         fprintf(stderr, "Note: Old option for generalized born electrostatics given:\n"
1235                 "Changing coulombtype from \"generalized-born\" to \"cut-off\" and instead\n"
1236                 "setting implicit-solvent value to \"GBSA\" in input section.\n");
1237     }
1238
1239     if (ir->sa_algorithm == esaSTILL)
1240     {
1241         sprintf(err_buf, "Still SA algorithm not available yet, use %s or %s instead\n", esa_names[esaAPPROX], esa_names[esaNO]);
1242         CHECK(ir->sa_algorithm == esaSTILL);
1243     }
1244
1245     if (ir->implicit_solvent == eisGBSA)
1246     {
1247         sprintf(err_buf, "With GBSA implicit solvent, rgbradii must be equal to rlist.");
1248         CHECK(ir->rgbradii != ir->rlist);
1249
1250         if (ir->coulombtype != eelCUT)
1251         {
1252             sprintf(err_buf, "With GBSA, coulombtype must be equal to %s\n", eel_names[eelCUT]);
1253             CHECK(ir->coulombtype != eelCUT);
1254         }
1255         if (ir->vdwtype != evdwCUT)
1256         {
1257             sprintf(err_buf, "With GBSA, vdw-type must be equal to %s\n", evdw_names[evdwCUT]);
1258             CHECK(ir->vdwtype != evdwCUT);
1259         }
1260         if (ir->nstgbradii < 1)
1261         {
1262             sprintf(warn_buf, "Using GBSA with nstgbradii<1, setting nstgbradii=1");
1263             warning_note(wi, warn_buf);
1264             ir->nstgbradii = 1;
1265         }
1266         if (ir->sa_algorithm == esaNO)
1267         {
1268             sprintf(warn_buf, "No SA (non-polar) calculation requested together with GB. Are you sure this is what you want?\n");
1269             warning_note(wi, warn_buf);
1270         }
1271         if (ir->sa_surface_tension < 0 && ir->sa_algorithm != esaNO)
1272         {
1273             sprintf(warn_buf, "Value of sa_surface_tension is < 0. Changing it to 2.05016 or 2.25936 kJ/nm^2/mol for Still and HCT/OBC respectively\n");
1274             warning_note(wi, warn_buf);
1275
1276             if (ir->gb_algorithm == egbSTILL)
1277             {
1278                 ir->sa_surface_tension = 0.0049 * CAL2JOULE * 100;
1279             }
1280             else
1281             {
1282                 ir->sa_surface_tension = 0.0054 * CAL2JOULE * 100;
1283             }
1284         }
1285         if (ir->sa_surface_tension == 0 && ir->sa_algorithm != esaNO)
1286         {
1287             sprintf(err_buf, "Surface tension set to 0 while SA-calculation requested\n");
1288             CHECK(ir->sa_surface_tension == 0 && ir->sa_algorithm != esaNO);
1289         }
1290
1291     }
1292
1293     if (ir->bAdress)
1294     {
1295         if (ir->cutoff_scheme != ecutsGROUP)
1296         {
1297             warning_error(wi, "AdresS simulation supports only cutoff-scheme=group");
1298         }
1299         if (!EI_SD(ir->eI))
1300         {
1301             warning_error(wi, "AdresS simulation supports only stochastic dynamics");
1302         }
1303         if (ir->epc != epcNO)
1304         {
1305             warning_error(wi, "AdresS simulation does not support pressure coupling");
1306         }
1307         if (EEL_FULL(ir->coulombtype))
1308         {
1309             warning_error(wi, "AdresS simulation does not support long-range electrostatics");
1310         }
1311     }
1312 }
1313
1314 /* count the number of text elemets separated by whitespace in a string.
1315     str = the input string
1316     maxptr = the maximum number of allowed elements
1317     ptr = the output array of pointers to the first character of each element
1318     returns: the number of elements. */
1319 int str_nelem(const char *str, int maxptr, char *ptr[])
1320 {
1321     int   np = 0;
1322     char *copy0, *copy;
1323
1324     copy0 = strdup(str);
1325     copy  = copy0;
1326     ltrim(copy);
1327     while (*copy != '\0')
1328     {
1329         if (np >= maxptr)
1330         {
1331             gmx_fatal(FARGS, "Too many groups on line: '%s' (max is %d)",
1332                       str, maxptr);
1333         }
1334         if (ptr)
1335         {
1336             ptr[np] = copy;
1337         }
1338         np++;
1339         while ((*copy != '\0') && !isspace(*copy))
1340         {
1341             copy++;
1342         }
1343         if (*copy != '\0')
1344         {
1345             *copy = '\0';
1346             copy++;
1347         }
1348         ltrim(copy);
1349     }
1350     if (ptr == NULL)
1351     {
1352         sfree(copy0);
1353     }
1354
1355     return np;
1356 }
1357
1358 /* interpret a number of doubles from a string and put them in an array,
1359    after allocating space for them.
1360    str = the input string
1361    n = the (pre-allocated) number of doubles read
1362    r = the output array of doubles. */
1363 static void parse_n_real(char *str, int *n, real **r)
1364 {
1365     char *ptr[MAXPTR];
1366     int   i;
1367
1368     *n = str_nelem(str, MAXPTR, ptr);
1369
1370     snew(*r, *n);
1371     for (i = 0; i < *n; i++)
1372     {
1373         (*r)[i] = strtod(ptr[i], NULL);
1374     }
1375 }
1376
1377 static void do_fep_params(t_inputrec *ir, char fep_lambda[][STRLEN], char weights[STRLEN])
1378 {
1379
1380     int         i, j, max_n_lambda, nweights, nfep[efptNR];
1381     t_lambda   *fep    = ir->fepvals;
1382     t_expanded *expand = ir->expandedvals;
1383     real      **count_fep_lambdas;
1384     gmx_bool    bOneLambda = TRUE;
1385
1386     snew(count_fep_lambdas, efptNR);
1387
1388     /* FEP input processing */
1389     /* first, identify the number of lambda values for each type.
1390        All that are nonzero must have the same number */
1391
1392     for (i = 0; i < efptNR; i++)
1393     {
1394         parse_n_real(fep_lambda[i], &(nfep[i]), &(count_fep_lambdas[i]));
1395     }
1396
1397     /* now, determine the number of components.  All must be either zero, or equal. */
1398
1399     max_n_lambda = 0;
1400     for (i = 0; i < efptNR; i++)
1401     {
1402         if (nfep[i] > max_n_lambda)
1403         {
1404             max_n_lambda = nfep[i];  /* here's a nonzero one.  All of them
1405                                         must have the same number if its not zero.*/
1406             break;
1407         }
1408     }
1409
1410     for (i = 0; i < efptNR; i++)
1411     {
1412         if (nfep[i] == 0)
1413         {
1414             ir->fepvals->separate_dvdl[i] = FALSE;
1415         }
1416         else if (nfep[i] == max_n_lambda)
1417         {
1418             if (i != efptTEMPERATURE)  /* we treat this differently -- not really a reason to compute the derivative with
1419                                           respect to the temperature currently */
1420             {
1421                 ir->fepvals->separate_dvdl[i] = TRUE;
1422             }
1423         }
1424         else
1425         {
1426             gmx_fatal(FARGS, "Number of lambdas (%d) for FEP type %s not equal to number of other types (%d)",
1427                       nfep[i], efpt_names[i], max_n_lambda);
1428         }
1429     }
1430     /* we don't print out dhdl if the temperature is changing, since we can't correctly define dhdl in this case */
1431     ir->fepvals->separate_dvdl[efptTEMPERATURE] = FALSE;
1432
1433     /* the number of lambdas is the number we've read in, which is either zero
1434        or the same for all */
1435     fep->n_lambda = max_n_lambda;
1436
1437     /* allocate space for the array of lambda values */
1438     snew(fep->all_lambda, efptNR);
1439     /* if init_lambda is defined, we need to set lambda */
1440     if ((fep->init_lambda > 0) && (fep->n_lambda == 0))
1441     {
1442         ir->fepvals->separate_dvdl[efptFEP] = TRUE;
1443     }
1444     /* otherwise allocate the space for all of the lambdas, and transfer the data */
1445     for (i = 0; i < efptNR; i++)
1446     {
1447         snew(fep->all_lambda[i], fep->n_lambda);
1448         if (nfep[i] > 0)  /* if it's zero, then the count_fep_lambda arrays
1449                              are zero */
1450         {
1451             for (j = 0; j < fep->n_lambda; j++)
1452             {
1453                 fep->all_lambda[i][j] = (double)count_fep_lambdas[i][j];
1454             }
1455             sfree(count_fep_lambdas[i]);
1456         }
1457     }
1458     sfree(count_fep_lambdas);
1459
1460     /* "fep-vals" is either zero or the full number. If zero, we'll need to define fep-lambdas for internal
1461        bookkeeping -- for now, init_lambda */
1462
1463     if ((nfep[efptFEP] == 0) && (fep->init_lambda >= 0))
1464     {
1465         for (i = 0; i < fep->n_lambda; i++)
1466         {
1467             fep->all_lambda[efptFEP][i] = fep->init_lambda;
1468         }
1469     }
1470
1471     /* check to see if only a single component lambda is defined, and soft core is defined.
1472        In this case, turn on coulomb soft core */
1473
1474     if (max_n_lambda == 0)
1475     {
1476         bOneLambda = TRUE;
1477     }
1478     else
1479     {
1480         for (i = 0; i < efptNR; i++)
1481         {
1482             if ((nfep[i] != 0) && (i != efptFEP))
1483             {
1484                 bOneLambda = FALSE;
1485             }
1486         }
1487     }
1488     if ((bOneLambda) && (fep->sc_alpha > 0))
1489     {
1490         fep->bScCoul = TRUE;
1491     }
1492
1493     /* Fill in the others with the efptFEP if they are not explicitly
1494        specified (i.e. nfep[i] == 0).  This means if fep is not defined,
1495        they are all zero. */
1496
1497     for (i = 0; i < efptNR; i++)
1498     {
1499         if ((nfep[i] == 0) && (i != efptFEP))
1500         {
1501             for (j = 0; j < fep->n_lambda; j++)
1502             {
1503                 fep->all_lambda[i][j] = fep->all_lambda[efptFEP][j];
1504             }
1505         }
1506     }
1507
1508
1509     /* make it easier if sc_r_power = 48 by increasing it to the 4th power, to be in the right scale. */
1510     if (fep->sc_r_power == 48)
1511     {
1512         if (fep->sc_alpha > 0.1)
1513         {
1514             gmx_fatal(FARGS, "sc_alpha (%f) for sc_r_power = 48 should usually be between 0.001 and 0.004", fep->sc_alpha);
1515         }
1516     }
1517
1518     expand = ir->expandedvals;
1519     /* now read in the weights */
1520     parse_n_real(weights, &nweights, &(expand->init_lambda_weights));
1521     if (nweights == 0)
1522     {
1523         expand->bInit_weights = FALSE;
1524         snew(expand->init_lambda_weights, fep->n_lambda); /* initialize to zero */
1525     }
1526     else if (nweights != fep->n_lambda)
1527     {
1528         gmx_fatal(FARGS, "Number of weights (%d) is not equal to number of lambda values (%d)",
1529                   nweights, fep->n_lambda);
1530     }
1531     else
1532     {
1533         expand->bInit_weights = TRUE;
1534     }
1535     if ((expand->nstexpanded < 0) && (ir->efep != efepNO))
1536     {
1537         expand->nstexpanded = fep->nstdhdl;
1538         /* if you don't specify nstexpanded when doing expanded ensemble free energy calcs, it is set to nstdhdl */
1539     }
1540     if ((expand->nstexpanded < 0) && ir->bSimTemp)
1541     {
1542         expand->nstexpanded = 2*(int)(ir->opts.tau_t[0]/ir->delta_t);
1543         /* if you don't specify nstexpanded when doing expanded ensemble simulated tempering, it is set to
1544            2*tau_t just to be careful so it's not to frequent  */
1545     }
1546 }
1547
1548
1549 static void do_simtemp_params(t_inputrec *ir)
1550 {
1551
1552     snew(ir->simtempvals->temperatures, ir->fepvals->n_lambda);
1553     GetSimTemps(ir->fepvals->n_lambda, ir->simtempvals, ir->fepvals->all_lambda[efptTEMPERATURE]);
1554
1555     return;
1556 }
1557
1558 static void do_wall_params(t_inputrec *ir,
1559                            char *wall_atomtype, char *wall_density,
1560                            t_gromppopts *opts)
1561 {
1562     int    nstr, i;
1563     char  *names[MAXPTR];
1564     double dbl;
1565
1566     opts->wall_atomtype[0] = NULL;
1567     opts->wall_atomtype[1] = NULL;
1568
1569     ir->wall_atomtype[0] = -1;
1570     ir->wall_atomtype[1] = -1;
1571     ir->wall_density[0]  = 0;
1572     ir->wall_density[1]  = 0;
1573
1574     if (ir->nwall > 0)
1575     {
1576         nstr = str_nelem(wall_atomtype, MAXPTR, names);
1577         if (nstr != ir->nwall)
1578         {
1579             gmx_fatal(FARGS, "Expected %d elements for wall_atomtype, found %d",
1580                       ir->nwall, nstr);
1581         }
1582         for (i = 0; i < ir->nwall; i++)
1583         {
1584             opts->wall_atomtype[i] = strdup(names[i]);
1585         }
1586
1587         if (ir->wall_type == ewt93 || ir->wall_type == ewt104)
1588         {
1589             nstr = str_nelem(wall_density, MAXPTR, names);
1590             if (nstr != ir->nwall)
1591             {
1592                 gmx_fatal(FARGS, "Expected %d elements for wall-density, found %d", ir->nwall, nstr);
1593             }
1594             for (i = 0; i < ir->nwall; i++)
1595             {
1596                 sscanf(names[i], "%lf", &dbl);
1597                 if (dbl <= 0)
1598                 {
1599                     gmx_fatal(FARGS, "wall-density[%d] = %f\n", i, dbl);
1600                 }
1601                 ir->wall_density[i] = dbl;
1602             }
1603         }
1604     }
1605 }
1606
1607 static void add_wall_energrps(gmx_groups_t *groups, int nwall, t_symtab *symtab)
1608 {
1609     int     i;
1610     t_grps *grps;
1611     char    str[STRLEN];
1612
1613     if (nwall > 0)
1614     {
1615         srenew(groups->grpname, groups->ngrpname+nwall);
1616         grps = &(groups->grps[egcENER]);
1617         srenew(grps->nm_ind, grps->nr+nwall);
1618         for (i = 0; i < nwall; i++)
1619         {
1620             sprintf(str, "wall%d", i);
1621             groups->grpname[groups->ngrpname] = put_symtab(symtab, str);
1622             grps->nm_ind[grps->nr++]          = groups->ngrpname++;
1623         }
1624     }
1625 }
1626
1627 void read_expandedparams(int *ninp_p, t_inpfile **inp_p,
1628                          t_expanded *expand, warninp_t wi)
1629 {
1630     int        ninp, nerror = 0;
1631     t_inpfile *inp;
1632
1633     ninp   = *ninp_p;
1634     inp    = *inp_p;
1635
1636     /* read expanded ensemble parameters */
1637     CCTYPE ("expanded ensemble variables");
1638     ITYPE ("nstexpanded", expand->nstexpanded, -1);
1639     EETYPE("lmc-stats", expand->elamstats, elamstats_names);
1640     EETYPE("lmc-move", expand->elmcmove, elmcmove_names);
1641     EETYPE("lmc-weights-equil", expand->elmceq, elmceq_names);
1642     ITYPE ("weight-equil-number-all-lambda", expand->equil_n_at_lam, -1);
1643     ITYPE ("weight-equil-number-samples", expand->equil_samples, -1);
1644     ITYPE ("weight-equil-number-steps", expand->equil_steps, -1);
1645     RTYPE ("weight-equil-wl-delta", expand->equil_wl_delta, -1);
1646     RTYPE ("weight-equil-count-ratio", expand->equil_ratio, -1);
1647     CCTYPE("Seed for Monte Carlo in lambda space");
1648     ITYPE ("lmc-seed", expand->lmc_seed, -1);
1649     RTYPE ("mc-temperature", expand->mc_temp, -1);
1650     ITYPE ("lmc-repeats", expand->lmc_repeats, 1);
1651     ITYPE ("lmc-gibbsdelta", expand->gibbsdeltalam, -1);
1652     ITYPE ("lmc-forced-nstart", expand->lmc_forced_nstart, 0);
1653     EETYPE("symmetrized-transition-matrix", expand->bSymmetrizedTMatrix, yesno_names);
1654     ITYPE("nst-transition-matrix", expand->nstTij, -1);
1655     ITYPE ("mininum-var-min", expand->minvarmin, 100); /*default is reasonable */
1656     ITYPE ("weight-c-range", expand->c_range, 0);      /* default is just C=0 */
1657     RTYPE ("wl-scale", expand->wl_scale, 0.8);
1658     RTYPE ("wl-ratio", expand->wl_ratio, 0.8);
1659     RTYPE ("init-wl-delta", expand->init_wl_delta, 1.0);
1660     EETYPE("wl-oneovert", expand->bWLoneovert, yesno_names);
1661
1662     *ninp_p   = ninp;
1663     *inp_p    = inp;
1664
1665     return;
1666 }
1667
1668 void get_ir(const char *mdparin, const char *mdparout,
1669             t_inputrec *ir, t_gromppopts *opts,
1670             warninp_t wi)
1671 {
1672     char       *dumstr[2];
1673     double      dumdub[2][6];
1674     t_inpfile  *inp;
1675     const char *tmp;
1676     int         i, j, m, ninp;
1677     char        warn_buf[STRLEN];
1678     t_lambda   *fep    = ir->fepvals;
1679     t_expanded *expand = ir->expandedvals;
1680
1681     inp = read_inpfile(mdparin, &ninp, wi);
1682
1683     snew(dumstr[0], STRLEN);
1684     snew(dumstr[1], STRLEN);
1685
1686     /* remove the following deprecated commands */
1687     REM_TYPE("title");
1688     REM_TYPE("cpp");
1689     REM_TYPE("domain-decomposition");
1690     REM_TYPE("andersen-seed");
1691     REM_TYPE("dihre");
1692     REM_TYPE("dihre-fc");
1693     REM_TYPE("dihre-tau");
1694     REM_TYPE("nstdihreout");
1695     REM_TYPE("nstcheckpoint");
1696
1697     /* replace the following commands with the clearer new versions*/
1698     REPL_TYPE("unconstrained-start", "continuation");
1699     REPL_TYPE("foreign-lambda", "fep-lambdas");
1700
1701     CCTYPE ("VARIOUS PREPROCESSING OPTIONS");
1702     CTYPE ("Preprocessor information: use cpp syntax.");
1703     CTYPE ("e.g.: -I/home/joe/doe -I/home/mary/roe");
1704     STYPE ("include", opts->include,  NULL);
1705     CTYPE ("e.g.: -DPOSRES -DFLEXIBLE (note these variable names are case sensitive)");
1706     STYPE ("define",  opts->define,   NULL);
1707
1708     CCTYPE ("RUN CONTROL PARAMETERS");
1709     EETYPE("integrator",  ir->eI,         ei_names);
1710     CTYPE ("Start time and timestep in ps");
1711     RTYPE ("tinit",   ir->init_t, 0.0);
1712     RTYPE ("dt",      ir->delta_t,    0.001);
1713     STEPTYPE ("nsteps",   ir->nsteps,     0);
1714     CTYPE ("For exact run continuation or redoing part of a run");
1715     STEPTYPE ("init-step", ir->init_step,  0);
1716     CTYPE ("Part index is updated automatically on checkpointing (keeps files separate)");
1717     ITYPE ("simulation-part", ir->simulation_part, 1);
1718     CTYPE ("mode for center of mass motion removal");
1719     EETYPE("comm-mode",   ir->comm_mode,  ecm_names);
1720     CTYPE ("number of steps for center of mass motion removal");
1721     ITYPE ("nstcomm", ir->nstcomm,    100);
1722     CTYPE ("group(s) for center of mass motion removal");
1723     STYPE ("comm-grps",   vcm,            NULL);
1724
1725     CCTYPE ("LANGEVIN DYNAMICS OPTIONS");
1726     CTYPE ("Friction coefficient (amu/ps) and random seed");
1727     RTYPE ("bd-fric",     ir->bd_fric,    0.0);
1728     ITYPE ("ld-seed",     ir->ld_seed,    1993);
1729
1730     /* Em stuff */
1731     CCTYPE ("ENERGY MINIMIZATION OPTIONS");
1732     CTYPE ("Force tolerance and initial step-size");
1733     RTYPE ("emtol",       ir->em_tol,     10.0);
1734     RTYPE ("emstep",      ir->em_stepsize, 0.01);
1735     CTYPE ("Max number of iterations in relax-shells");
1736     ITYPE ("niter",       ir->niter,      20);
1737     CTYPE ("Step size (ps^2) for minimization of flexible constraints");
1738     RTYPE ("fcstep",      ir->fc_stepsize, 0);
1739     CTYPE ("Frequency of steepest descents steps when doing CG");
1740     ITYPE ("nstcgsteep",  ir->nstcgsteep, 1000);
1741     ITYPE ("nbfgscorr",   ir->nbfgscorr,  10);
1742
1743     CCTYPE ("TEST PARTICLE INSERTION OPTIONS");
1744     RTYPE ("rtpi",    ir->rtpi,   0.05);
1745
1746     /* Output options */
1747     CCTYPE ("OUTPUT CONTROL OPTIONS");
1748     CTYPE ("Output frequency for coords (x), velocities (v) and forces (f)");
1749     ITYPE ("nstxout", ir->nstxout,    0);
1750     ITYPE ("nstvout", ir->nstvout,    0);
1751     ITYPE ("nstfout", ir->nstfout,    0);
1752     ir->nstcheckpoint = 1000;
1753     CTYPE ("Output frequency for energies to log file and energy file");
1754     ITYPE ("nstlog",  ir->nstlog, 1000);
1755     ITYPE ("nstcalcenergy", ir->nstcalcenergy, 100);
1756     ITYPE ("nstenergy",   ir->nstenergy,  1000);
1757     CTYPE ("Output frequency and precision for .xtc file");
1758     ITYPE ("nstxtcout",   ir->nstxtcout,  0);
1759     RTYPE ("xtc-precision", ir->xtcprec,   1000.0);
1760     CTYPE ("This selects the subset of atoms for the .xtc file. You can");
1761     CTYPE ("select multiple groups. By default all atoms will be written.");
1762     STYPE ("xtc-grps",    xtc_grps,       NULL);
1763     CTYPE ("Selection of energy groups");
1764     STYPE ("energygrps",  energy,         NULL);
1765
1766     /* Neighbor searching */
1767     CCTYPE ("NEIGHBORSEARCHING PARAMETERS");
1768     CTYPE ("cut-off scheme (group: using charge groups, Verlet: particle based cut-offs)");
1769     EETYPE("cutoff-scheme",     ir->cutoff_scheme,    ecutscheme_names);
1770     CTYPE ("nblist update frequency");
1771     ITYPE ("nstlist", ir->nstlist,    10);
1772     CTYPE ("ns algorithm (simple or grid)");
1773     EETYPE("ns-type",     ir->ns_type,    ens_names);
1774     /* set ndelta to the optimal value of 2 */
1775     ir->ndelta = 2;
1776     CTYPE ("Periodic boundary conditions: xyz, no, xy");
1777     EETYPE("pbc",         ir->ePBC,       epbc_names);
1778     EETYPE("periodic-molecules", ir->bPeriodicMols, yesno_names);
1779     CTYPE ("Allowed energy drift due to the Verlet buffer in kJ/mol/ps per atom,");
1780     CTYPE ("a value of -1 means: use rlist");
1781     RTYPE("verlet-buffer-drift", ir->verletbuf_drift,    0.005);
1782     CTYPE ("nblist cut-off");
1783     RTYPE ("rlist",   ir->rlist,  1.0);
1784     CTYPE ("long-range cut-off for switched potentials");
1785     RTYPE ("rlistlong",   ir->rlistlong,  -1);
1786     ITYPE ("nstcalclr",   ir->nstcalclr,  -1);
1787
1788     /* Electrostatics */
1789     CCTYPE ("OPTIONS FOR ELECTROSTATICS AND VDW");
1790     CTYPE ("Method for doing electrostatics");
1791     EETYPE("coulombtype", ir->coulombtype,    eel_names);
1792     EETYPE("coulomb-modifier",    ir->coulomb_modifier,    eintmod_names);
1793     CTYPE ("cut-off lengths");
1794     RTYPE ("rcoulomb-switch", ir->rcoulomb_switch,    0.0);
1795     RTYPE ("rcoulomb",    ir->rcoulomb,   1.0);
1796     CTYPE ("Relative dielectric constant for the medium and the reaction field");
1797     RTYPE ("epsilon-r",   ir->epsilon_r,  1.0);
1798     RTYPE ("epsilon-rf",  ir->epsilon_rf, 0.0);
1799     CTYPE ("Method for doing Van der Waals");
1800     EETYPE("vdw-type",    ir->vdwtype,    evdw_names);
1801     EETYPE("vdw-modifier",    ir->vdw_modifier,    eintmod_names);
1802     CTYPE ("cut-off lengths");
1803     RTYPE ("rvdw-switch", ir->rvdw_switch,    0.0);
1804     RTYPE ("rvdw",    ir->rvdw,   1.0);
1805     CTYPE ("Apply long range dispersion corrections for Energy and Pressure");
1806     EETYPE("DispCorr",    ir->eDispCorr,  edispc_names);
1807     CTYPE ("Extension of the potential lookup tables beyond the cut-off");
1808     RTYPE ("table-extension", ir->tabext, 1.0);
1809     CTYPE ("Separate tables between energy group pairs");
1810     STYPE ("energygrp-table", egptable,   NULL);
1811     CTYPE ("Spacing for the PME/PPPM FFT grid");
1812     RTYPE ("fourierspacing", ir->fourier_spacing, 0.12);
1813     CTYPE ("FFT grid size, when a value is 0 fourierspacing will be used");
1814     ITYPE ("fourier-nx",  ir->nkx,         0);
1815     ITYPE ("fourier-ny",  ir->nky,         0);
1816     ITYPE ("fourier-nz",  ir->nkz,         0);
1817     CTYPE ("EWALD/PME/PPPM parameters");
1818     ITYPE ("pme-order",   ir->pme_order,   4);
1819     RTYPE ("ewald-rtol",  ir->ewald_rtol, 0.00001);
1820     EETYPE("ewald-geometry", ir->ewald_geometry, eewg_names);
1821     RTYPE ("epsilon-surface", ir->epsilon_surface, 0.0);
1822     EETYPE("optimize-fft", ir->bOptFFT,  yesno_names);
1823
1824     CCTYPE("IMPLICIT SOLVENT ALGORITHM");
1825     EETYPE("implicit-solvent", ir->implicit_solvent, eis_names);
1826
1827     CCTYPE ("GENERALIZED BORN ELECTROSTATICS");
1828     CTYPE ("Algorithm for calculating Born radii");
1829     EETYPE("gb-algorithm", ir->gb_algorithm, egb_names);
1830     CTYPE ("Frequency of calculating the Born radii inside rlist");
1831     ITYPE ("nstgbradii", ir->nstgbradii, 1);
1832     CTYPE ("Cutoff for Born radii calculation; the contribution from atoms");
1833     CTYPE ("between rlist and rgbradii is updated every nstlist steps");
1834     RTYPE ("rgbradii",  ir->rgbradii, 1.0);
1835     CTYPE ("Dielectric coefficient of the implicit solvent");
1836     RTYPE ("gb-epsilon-solvent", ir->gb_epsilon_solvent, 80.0);
1837     CTYPE ("Salt concentration in M for Generalized Born models");
1838     RTYPE ("gb-saltconc",  ir->gb_saltconc, 0.0);
1839     CTYPE ("Scaling factors used in the OBC GB model. Default values are OBC(II)");
1840     RTYPE ("gb-obc-alpha", ir->gb_obc_alpha, 1.0);
1841     RTYPE ("gb-obc-beta", ir->gb_obc_beta, 0.8);
1842     RTYPE ("gb-obc-gamma", ir->gb_obc_gamma, 4.85);
1843     RTYPE ("gb-dielectric-offset", ir->gb_dielectric_offset, 0.009);
1844     EETYPE("sa-algorithm", ir->sa_algorithm, esa_names);
1845     CTYPE ("Surface tension (kJ/mol/nm^2) for the SA (nonpolar surface) part of GBSA");
1846     CTYPE ("The value -1 will set default value for Still/HCT/OBC GB-models.");
1847     RTYPE ("sa-surface-tension", ir->sa_surface_tension, -1);
1848
1849     /* Coupling stuff */
1850     CCTYPE ("OPTIONS FOR WEAK COUPLING ALGORITHMS");
1851     CTYPE ("Temperature coupling");
1852     EETYPE("tcoupl",  ir->etc,        etcoupl_names);
1853     ITYPE ("nsttcouple", ir->nsttcouple,  -1);
1854     ITYPE("nh-chain-length",     ir->opts.nhchainlength, NHCHAINLENGTH);
1855     EETYPE("print-nose-hoover-chain-variables", ir->bPrintNHChains, yesno_names);
1856     CTYPE ("Groups to couple separately");
1857     STYPE ("tc-grps",     tcgrps,         NULL);
1858     CTYPE ("Time constant (ps) and reference temperature (K)");
1859     STYPE ("tau-t",   tau_t,      NULL);
1860     STYPE ("ref-t",   ref_t,      NULL);
1861     CTYPE ("pressure coupling");
1862     EETYPE("pcoupl",  ir->epc,        epcoupl_names);
1863     EETYPE("pcoupltype",  ir->epct,       epcoupltype_names);
1864     ITYPE ("nstpcouple", ir->nstpcouple,  -1);
1865     CTYPE ("Time constant (ps), compressibility (1/bar) and reference P (bar)");
1866     RTYPE ("tau-p",   ir->tau_p,  1.0);
1867     STYPE ("compressibility", dumstr[0],  NULL);
1868     STYPE ("ref-p",       dumstr[1],      NULL);
1869     CTYPE ("Scaling of reference coordinates, No, All or COM");
1870     EETYPE ("refcoord-scaling", ir->refcoord_scaling, erefscaling_names);
1871
1872     /* QMMM */
1873     CCTYPE ("OPTIONS FOR QMMM calculations");
1874     EETYPE("QMMM", ir->bQMMM, yesno_names);
1875     CTYPE ("Groups treated Quantum Mechanically");
1876     STYPE ("QMMM-grps",  QMMM,          NULL);
1877     CTYPE ("QM method");
1878     STYPE("QMmethod",     QMmethod, NULL);
1879     CTYPE ("QMMM scheme");
1880     EETYPE("QMMMscheme",  ir->QMMMscheme,    eQMMMscheme_names);
1881     CTYPE ("QM basisset");
1882     STYPE("QMbasis",      QMbasis, NULL);
1883     CTYPE ("QM charge");
1884     STYPE ("QMcharge",    QMcharge, NULL);
1885     CTYPE ("QM multiplicity");
1886     STYPE ("QMmult",      QMmult, NULL);
1887     CTYPE ("Surface Hopping");
1888     STYPE ("SH",          bSH, NULL);
1889     CTYPE ("CAS space options");
1890     STYPE ("CASorbitals",      CASorbitals,   NULL);
1891     STYPE ("CASelectrons",     CASelectrons,  NULL);
1892     STYPE ("SAon", SAon, NULL);
1893     STYPE ("SAoff", SAoff, NULL);
1894     STYPE ("SAsteps",  SAsteps, NULL);
1895     CTYPE ("Scale factor for MM charges");
1896     RTYPE ("MMChargeScaleFactor", ir->scalefactor, 1.0);
1897     CTYPE ("Optimization of QM subsystem");
1898     STYPE ("bOPT",          bOPT, NULL);
1899     STYPE ("bTS",          bTS, NULL);
1900
1901     /* Simulated annealing */
1902     CCTYPE("SIMULATED ANNEALING");
1903     CTYPE ("Type of annealing for each temperature group (no/single/periodic)");
1904     STYPE ("annealing",   anneal,      NULL);
1905     CTYPE ("Number of time points to use for specifying annealing in each group");
1906     STYPE ("annealing-npoints", anneal_npoints, NULL);
1907     CTYPE ("List of times at the annealing points for each group");
1908     STYPE ("annealing-time",       anneal_time,       NULL);
1909     CTYPE ("Temp. at each annealing point, for each group.");
1910     STYPE ("annealing-temp",  anneal_temp,  NULL);
1911
1912     /* Startup run */
1913     CCTYPE ("GENERATE VELOCITIES FOR STARTUP RUN");
1914     EETYPE("gen-vel",     opts->bGenVel,  yesno_names);
1915     RTYPE ("gen-temp",    opts->tempi,    300.0);
1916     ITYPE ("gen-seed",    opts->seed,     173529);
1917
1918     /* Shake stuff */
1919     CCTYPE ("OPTIONS FOR BONDS");
1920     EETYPE("constraints", opts->nshake,   constraints);
1921     CTYPE ("Type of constraint algorithm");
1922     EETYPE("constraint-algorithm",  ir->eConstrAlg, econstr_names);
1923     CTYPE ("Do not constrain the start configuration");
1924     EETYPE("continuation", ir->bContinuation, yesno_names);
1925     CTYPE ("Use successive overrelaxation to reduce the number of shake iterations");
1926     EETYPE("Shake-SOR", ir->bShakeSOR, yesno_names);
1927     CTYPE ("Relative tolerance of shake");
1928     RTYPE ("shake-tol", ir->shake_tol, 0.0001);
1929     CTYPE ("Highest order in the expansion of the constraint coupling matrix");
1930     ITYPE ("lincs-order", ir->nProjOrder, 4);
1931     CTYPE ("Number of iterations in the final step of LINCS. 1 is fine for");
1932     CTYPE ("normal simulations, but use 2 to conserve energy in NVE runs.");
1933     CTYPE ("For energy minimization with constraints it should be 4 to 8.");
1934     ITYPE ("lincs-iter", ir->nLincsIter, 1);
1935     CTYPE ("Lincs will write a warning to the stderr if in one step a bond");
1936     CTYPE ("rotates over more degrees than");
1937     RTYPE ("lincs-warnangle", ir->LincsWarnAngle, 30.0);
1938     CTYPE ("Convert harmonic bonds to morse potentials");
1939     EETYPE("morse",       opts->bMorse, yesno_names);
1940
1941     /* Energy group exclusions */
1942     CCTYPE ("ENERGY GROUP EXCLUSIONS");
1943     CTYPE ("Pairs of energy groups for which all non-bonded interactions are excluded");
1944     STYPE ("energygrp-excl", egpexcl,     NULL);
1945
1946     /* Walls */
1947     CCTYPE ("WALLS");
1948     CTYPE ("Number of walls, type, atom types, densities and box-z scale factor for Ewald");
1949     ITYPE ("nwall", ir->nwall, 0);
1950     EETYPE("wall-type",     ir->wall_type,   ewt_names);
1951     RTYPE ("wall-r-linpot", ir->wall_r_linpot, -1);
1952     STYPE ("wall-atomtype", wall_atomtype, NULL);
1953     STYPE ("wall-density",  wall_density,  NULL);
1954     RTYPE ("wall-ewald-zfac", ir->wall_ewald_zfac, 3);
1955
1956     /* COM pulling */
1957     CCTYPE("COM PULLING");
1958     CTYPE("Pull type: no, umbrella, constraint or constant-force");
1959     EETYPE("pull",          ir->ePull, epull_names);
1960     if (ir->ePull != epullNO)
1961     {
1962         snew(ir->pull, 1);
1963         pull_grp = read_pullparams(&ninp, &inp, ir->pull, &opts->pull_start, wi);
1964     }
1965
1966     /* Enforced rotation */
1967     CCTYPE("ENFORCED ROTATION");
1968     CTYPE("Enforced rotation: No or Yes");
1969     EETYPE("rotation",       ir->bRot, yesno_names);
1970     if (ir->bRot)
1971     {
1972         snew(ir->rot, 1);
1973         rot_grp = read_rotparams(&ninp, &inp, ir->rot, wi);
1974     }
1975
1976     /* Refinement */
1977     CCTYPE("NMR refinement stuff");
1978     CTYPE ("Distance restraints type: No, Simple or Ensemble");
1979     EETYPE("disre",       ir->eDisre,     edisre_names);
1980     CTYPE ("Force weighting of pairs in one distance restraint: Conservative or Equal");
1981     EETYPE("disre-weighting", ir->eDisreWeighting, edisreweighting_names);
1982     CTYPE ("Use sqrt of the time averaged times the instantaneous violation");
1983     EETYPE("disre-mixed", ir->bDisreMixed, yesno_names);
1984     RTYPE ("disre-fc",    ir->dr_fc,  1000.0);
1985     RTYPE ("disre-tau",   ir->dr_tau, 0.0);
1986     CTYPE ("Output frequency for pair distances to energy file");
1987     ITYPE ("nstdisreout", ir->nstdisreout, 100);
1988     CTYPE ("Orientation restraints: No or Yes");
1989     EETYPE("orire",       opts->bOrire,   yesno_names);
1990     CTYPE ("Orientation restraints force constant and tau for time averaging");
1991     RTYPE ("orire-fc",    ir->orires_fc,  0.0);
1992     RTYPE ("orire-tau",   ir->orires_tau, 0.0);
1993     STYPE ("orire-fitgrp", orirefitgrp,    NULL);
1994     CTYPE ("Output frequency for trace(SD) and S to energy file");
1995     ITYPE ("nstorireout", ir->nstorireout, 100);
1996
1997     /* free energy variables */
1998     CCTYPE ("Free energy variables");
1999     EETYPE("free-energy", ir->efep, efep_names);
2000     STYPE ("couple-moltype",  couple_moltype,  NULL);
2001     EETYPE("couple-lambda0", opts->couple_lam0, couple_lam);
2002     EETYPE("couple-lambda1", opts->couple_lam1, couple_lam);
2003     EETYPE("couple-intramol", opts->bCoupleIntra, yesno_names);
2004
2005     RTYPE ("init-lambda", fep->init_lambda, -1); /* start with -1 so
2006                                                     we can recognize if
2007                                                     it was not entered */
2008     ITYPE ("init-lambda-state", fep->init_fep_state, -1);
2009     RTYPE ("delta-lambda", fep->delta_lambda, 0.0);
2010     ITYPE ("nstdhdl", fep->nstdhdl, 50);
2011     STYPE ("fep-lambdas", fep_lambda[efptFEP], NULL);
2012     STYPE ("mass-lambdas", fep_lambda[efptMASS], NULL);
2013     STYPE ("coul-lambdas", fep_lambda[efptCOUL], NULL);
2014     STYPE ("vdw-lambdas", fep_lambda[efptVDW], NULL);
2015     STYPE ("bonded-lambdas", fep_lambda[efptBONDED], NULL);
2016     STYPE ("restraint-lambdas", fep_lambda[efptRESTRAINT], NULL);
2017     STYPE ("temperature-lambdas", fep_lambda[efptTEMPERATURE], NULL);
2018     ITYPE ("calc-lambda-neighbors", fep->lambda_neighbors, 1);
2019     STYPE ("init-lambda-weights", lambda_weights, NULL);
2020     EETYPE("dhdl-print-energy", fep->bPrintEnergy, yesno_names);
2021     RTYPE ("sc-alpha", fep->sc_alpha, 0.0);
2022     ITYPE ("sc-power", fep->sc_power, 1);
2023     RTYPE ("sc-r-power", fep->sc_r_power, 6.0);
2024     RTYPE ("sc-sigma", fep->sc_sigma, 0.3);
2025     EETYPE("sc-coul", fep->bScCoul, yesno_names);
2026     ITYPE ("dh_hist_size", fep->dh_hist_size, 0);
2027     RTYPE ("dh_hist_spacing", fep->dh_hist_spacing, 0.1);
2028     EETYPE("separate-dhdl-file", fep->separate_dhdl_file,
2029            separate_dhdl_file_names);
2030     EETYPE("dhdl-derivatives", fep->dhdl_derivatives, dhdl_derivatives_names);
2031     ITYPE ("dh_hist_size", fep->dh_hist_size, 0);
2032     RTYPE ("dh_hist_spacing", fep->dh_hist_spacing, 0.1);
2033
2034     /* Non-equilibrium MD stuff */
2035     CCTYPE("Non-equilibrium MD stuff");
2036     STYPE ("acc-grps",    accgrps,        NULL);
2037     STYPE ("accelerate",  acc,            NULL);
2038     STYPE ("freezegrps",  freeze,         NULL);
2039     STYPE ("freezedim",   frdim,          NULL);
2040     RTYPE ("cos-acceleration", ir->cos_accel, 0);
2041     STYPE ("deform",      deform,         NULL);
2042
2043     /* simulated tempering variables */
2044     CCTYPE("simulated tempering variables");
2045     EETYPE("simulated-tempering", ir->bSimTemp, yesno_names);
2046     EETYPE("simulated-tempering-scaling", ir->simtempvals->eSimTempScale, esimtemp_names);
2047     RTYPE("sim-temp-low", ir->simtempvals->simtemp_low, 300.0);
2048     RTYPE("sim-temp-high", ir->simtempvals->simtemp_high, 300.0);
2049
2050     /* expanded ensemble variables */
2051     if (ir->efep == efepEXPANDED || ir->bSimTemp)
2052     {
2053         read_expandedparams(&ninp, &inp, expand, wi);
2054     }
2055
2056     /* Electric fields */
2057     CCTYPE("Electric fields");
2058     CTYPE ("Format is number of terms (int) and for all terms an amplitude (real)");
2059     CTYPE ("and a phase angle (real)");
2060     STYPE ("E-x",     efield_x,   NULL);
2061     STYPE ("E-xt",    efield_xt,  NULL);
2062     STYPE ("E-y",     efield_y,   NULL);
2063     STYPE ("E-yt",    efield_yt,  NULL);
2064     STYPE ("E-z",     efield_z,   NULL);
2065     STYPE ("E-zt",    efield_zt,  NULL);
2066
2067     /* AdResS defined thingies */
2068     CCTYPE ("AdResS parameters");
2069     EETYPE("adress",       ir->bAdress, yesno_names);
2070     if (ir->bAdress)
2071     {
2072         snew(ir->adress, 1);
2073         read_adressparams(&ninp, &inp, ir->adress, wi);
2074     }
2075
2076     /* User defined thingies */
2077     CCTYPE ("User defined thingies");
2078     STYPE ("user1-grps",  user1,          NULL);
2079     STYPE ("user2-grps",  user2,          NULL);
2080     ITYPE ("userint1",    ir->userint1,   0);
2081     ITYPE ("userint2",    ir->userint2,   0);
2082     ITYPE ("userint3",    ir->userint3,   0);
2083     ITYPE ("userint4",    ir->userint4,   0);
2084     RTYPE ("userreal1",   ir->userreal1,  0);
2085     RTYPE ("userreal2",   ir->userreal2,  0);
2086     RTYPE ("userreal3",   ir->userreal3,  0);
2087     RTYPE ("userreal4",   ir->userreal4,  0);
2088 #undef CTYPE
2089
2090     write_inpfile(mdparout, ninp, inp, FALSE, wi);
2091     for (i = 0; (i < ninp); i++)
2092     {
2093         sfree(inp[i].name);
2094         sfree(inp[i].value);
2095     }
2096     sfree(inp);
2097
2098     /* Process options if necessary */
2099     for (m = 0; m < 2; m++)
2100     {
2101         for (i = 0; i < 2*DIM; i++)
2102         {
2103             dumdub[m][i] = 0.0;
2104         }
2105         if (ir->epc)
2106         {
2107             switch (ir->epct)
2108             {
2109                 case epctISOTROPIC:
2110                     if (sscanf(dumstr[m], "%lf", &(dumdub[m][XX])) != 1)
2111                     {
2112                         warning_error(wi, "Pressure coupling not enough values (I need 1)");
2113                     }
2114                     dumdub[m][YY] = dumdub[m][ZZ] = dumdub[m][XX];
2115                     break;
2116                 case epctSEMIISOTROPIC:
2117                 case epctSURFACETENSION:
2118                     if (sscanf(dumstr[m], "%lf%lf",
2119                                &(dumdub[m][XX]), &(dumdub[m][ZZ])) != 2)
2120                     {
2121                         warning_error(wi, "Pressure coupling not enough values (I need 2)");
2122                     }
2123                     dumdub[m][YY] = dumdub[m][XX];
2124                     break;
2125                 case epctANISOTROPIC:
2126                     if (sscanf(dumstr[m], "%lf%lf%lf%lf%lf%lf",
2127                                &(dumdub[m][XX]), &(dumdub[m][YY]), &(dumdub[m][ZZ]),
2128                                &(dumdub[m][3]), &(dumdub[m][4]), &(dumdub[m][5])) != 6)
2129                     {
2130                         warning_error(wi, "Pressure coupling not enough values (I need 6)");
2131                     }
2132                     break;
2133                 default:
2134                     gmx_fatal(FARGS, "Pressure coupling type %s not implemented yet",
2135                               epcoupltype_names[ir->epct]);
2136             }
2137         }
2138     }
2139     clear_mat(ir->ref_p);
2140     clear_mat(ir->compress);
2141     for (i = 0; i < DIM; i++)
2142     {
2143         ir->ref_p[i][i]    = dumdub[1][i];
2144         ir->compress[i][i] = dumdub[0][i];
2145     }
2146     if (ir->epct == epctANISOTROPIC)
2147     {
2148         ir->ref_p[XX][YY] = dumdub[1][3];
2149         ir->ref_p[XX][ZZ] = dumdub[1][4];
2150         ir->ref_p[YY][ZZ] = dumdub[1][5];
2151         if (ir->ref_p[XX][YY] != 0 && ir->ref_p[XX][ZZ] != 0 && ir->ref_p[YY][ZZ] != 0)
2152         {
2153             warning(wi, "All off-diagonal reference pressures are non-zero. Are you sure you want to apply a threefold shear stress?\n");
2154         }
2155         ir->compress[XX][YY] = dumdub[0][3];
2156         ir->compress[XX][ZZ] = dumdub[0][4];
2157         ir->compress[YY][ZZ] = dumdub[0][5];
2158         for (i = 0; i < DIM; i++)
2159         {
2160             for (m = 0; m < i; m++)
2161             {
2162                 ir->ref_p[i][m]    = ir->ref_p[m][i];
2163                 ir->compress[i][m] = ir->compress[m][i];
2164             }
2165         }
2166     }
2167
2168     if (ir->comm_mode == ecmNO)
2169     {
2170         ir->nstcomm = 0;
2171     }
2172
2173     opts->couple_moltype = NULL;
2174     if (strlen(couple_moltype) > 0)
2175     {
2176         if (ir->efep != efepNO)
2177         {
2178             opts->couple_moltype = strdup(couple_moltype);
2179             if (opts->couple_lam0 == opts->couple_lam1)
2180             {
2181                 warning(wi, "The lambda=0 and lambda=1 states for coupling are identical");
2182             }
2183             if (ir->eI == eiMD && (opts->couple_lam0 == ecouplamNONE ||
2184                                    opts->couple_lam1 == ecouplamNONE))
2185             {
2186                 warning(wi, "For proper sampling of the (nearly) decoupled state, stochastic dynamics should be used");
2187             }
2188         }
2189         else
2190         {
2191             warning(wi, "Can not couple a molecule with free_energy = no");
2192         }
2193     }
2194     /* FREE ENERGY AND EXPANDED ENSEMBLE OPTIONS */
2195     if (ir->efep != efepNO)
2196     {
2197         if (fep->delta_lambda > 0)
2198         {
2199             ir->efep = efepSLOWGROWTH;
2200         }
2201     }
2202
2203     if (ir->bSimTemp)
2204     {
2205         fep->bPrintEnergy = TRUE;
2206         /* always print out the energy to dhdl if we are doing expanded ensemble, since we need the total energy
2207            if the temperature is changing. */
2208     }
2209
2210     if ((ir->efep != efepNO) || ir->bSimTemp)
2211     {
2212         ir->bExpanded = FALSE;
2213         if ((ir->efep == efepEXPANDED) || ir->bSimTemp)
2214         {
2215             ir->bExpanded = TRUE;
2216         }
2217         do_fep_params(ir, fep_lambda, lambda_weights);
2218         if (ir->bSimTemp) /* done after fep params */
2219         {
2220             do_simtemp_params(ir);
2221         }
2222     }
2223     else
2224     {
2225         ir->fepvals->n_lambda = 0;
2226     }
2227
2228     /* WALL PARAMETERS */
2229
2230     do_wall_params(ir, wall_atomtype, wall_density, opts);
2231
2232     /* ORIENTATION RESTRAINT PARAMETERS */
2233
2234     if (opts->bOrire && str_nelem(orirefitgrp, MAXPTR, NULL) != 1)
2235     {
2236         warning_error(wi, "ERROR: Need one orientation restraint fit group\n");
2237     }
2238
2239     /* DEFORMATION PARAMETERS */
2240
2241     clear_mat(ir->deform);
2242     for (i = 0; i < 6; i++)
2243     {
2244         dumdub[0][i] = 0;
2245     }
2246     m = sscanf(deform, "%lf %lf %lf %lf %lf %lf",
2247                &(dumdub[0][0]), &(dumdub[0][1]), &(dumdub[0][2]),
2248                &(dumdub[0][3]), &(dumdub[0][4]), &(dumdub[0][5]));
2249     for (i = 0; i < 3; i++)
2250     {
2251         ir->deform[i][i] = dumdub[0][i];
2252     }
2253     ir->deform[YY][XX] = dumdub[0][3];
2254     ir->deform[ZZ][XX] = dumdub[0][4];
2255     ir->deform[ZZ][YY] = dumdub[0][5];
2256     if (ir->epc != epcNO)
2257     {
2258         for (i = 0; i < 3; i++)
2259         {
2260             for (j = 0; j <= i; j++)
2261             {
2262                 if (ir->deform[i][j] != 0 && ir->compress[i][j] != 0)
2263                 {
2264                     warning_error(wi, "A box element has deform set and compressibility > 0");
2265                 }
2266             }
2267         }
2268         for (i = 0; i < 3; i++)
2269         {
2270             for (j = 0; j < i; j++)
2271             {
2272                 if (ir->deform[i][j] != 0)
2273                 {
2274                     for (m = j; m < DIM; m++)
2275                     {
2276                         if (ir->compress[m][j] != 0)
2277                         {
2278                             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.");
2279                             warning(wi, warn_buf);
2280                         }
2281                     }
2282                 }
2283             }
2284         }
2285     }
2286
2287     sfree(dumstr[0]);
2288     sfree(dumstr[1]);
2289 }
2290
2291 static int search_QMstring(char *s, int ng, const char *gn[])
2292 {
2293     /* same as normal search_string, but this one searches QM strings */
2294     int i;
2295
2296     for (i = 0; (i < ng); i++)
2297     {
2298         if (gmx_strcasecmp(s, gn[i]) == 0)
2299         {
2300             return i;
2301         }
2302     }
2303
2304     gmx_fatal(FARGS, "this QM method or basisset (%s) is not implemented\n!", s);
2305
2306     return -1;
2307
2308 } /* search_QMstring */
2309
2310
2311 int search_string(char *s, int ng, char *gn[])
2312 {
2313     int i;
2314
2315     for (i = 0; (i < ng); i++)
2316     {
2317         if (gmx_strcasecmp(s, gn[i]) == 0)
2318         {
2319             return i;
2320         }
2321     }
2322
2323     gmx_fatal(FARGS,
2324               "Group %s referenced in the .mdp file was not found in the index file.\n"
2325               "Group names must match either [moleculetype] names or custom index group\n"
2326               "names, in which case you must supply an index file to the '-n' option\n"
2327               "of grompp.",
2328               s);
2329
2330     return -1;
2331 }
2332
2333 static gmx_bool do_numbering(int natoms, gmx_groups_t *groups, int ng, char *ptrs[],
2334                              t_blocka *block, char *gnames[],
2335                              int gtype, int restnm,
2336                              int grptp, gmx_bool bVerbose,
2337                              warninp_t wi)
2338 {
2339     unsigned short *cbuf;
2340     t_grps         *grps = &(groups->grps[gtype]);
2341     int             i, j, gid, aj, ognr, ntot = 0;
2342     const char     *title;
2343     gmx_bool        bRest;
2344     char            warn_buf[STRLEN];
2345
2346     if (debug)
2347     {
2348         fprintf(debug, "Starting numbering %d groups of type %d\n", ng, gtype);
2349     }
2350
2351     title = gtypes[gtype];
2352
2353     snew(cbuf, natoms);
2354     /* Mark all id's as not set */
2355     for (i = 0; (i < natoms); i++)
2356     {
2357         cbuf[i] = NOGID;
2358     }
2359
2360     snew(grps->nm_ind, ng+1); /* +1 for possible rest group */
2361     for (i = 0; (i < ng); i++)
2362     {
2363         /* Lookup the group name in the block structure */
2364         gid = search_string(ptrs[i], block->nr, gnames);
2365         if ((grptp != egrptpONE) || (i == 0))
2366         {
2367             grps->nm_ind[grps->nr++] = gid;
2368         }
2369         if (debug)
2370         {
2371             fprintf(debug, "Found gid %d for group %s\n", gid, ptrs[i]);
2372         }
2373
2374         /* Now go over the atoms in the group */
2375         for (j = block->index[gid]; (j < block->index[gid+1]); j++)
2376         {
2377
2378             aj = block->a[j];
2379
2380             /* Range checking */
2381             if ((aj < 0) || (aj >= natoms))
2382             {
2383                 gmx_fatal(FARGS, "Invalid atom number %d in indexfile", aj);
2384             }
2385             /* Lookup up the old group number */
2386             ognr = cbuf[aj];
2387             if (ognr != NOGID)
2388             {
2389                 gmx_fatal(FARGS, "Atom %d in multiple %s groups (%d and %d)",
2390                           aj+1, title, ognr+1, i+1);
2391             }
2392             else
2393             {
2394                 /* Store the group number in buffer */
2395                 if (grptp == egrptpONE)
2396                 {
2397                     cbuf[aj] = 0;
2398                 }
2399                 else
2400                 {
2401                     cbuf[aj] = i;
2402                 }
2403                 ntot++;
2404             }
2405         }
2406     }
2407
2408     /* Now check whether we have done all atoms */
2409     bRest = FALSE;
2410     if (ntot != natoms)
2411     {
2412         if (grptp == egrptpALL)
2413         {
2414             gmx_fatal(FARGS, "%d atoms are not part of any of the %s groups",
2415                       natoms-ntot, title);
2416         }
2417         else if (grptp == egrptpPART)
2418         {
2419             sprintf(warn_buf, "%d atoms are not part of any of the %s groups",
2420                     natoms-ntot, title);
2421             warning_note(wi, warn_buf);
2422         }
2423         /* Assign all atoms currently unassigned to a rest group */
2424         for (j = 0; (j < natoms); j++)
2425         {
2426             if (cbuf[j] == NOGID)
2427             {
2428                 cbuf[j] = grps->nr;
2429                 bRest   = TRUE;
2430             }
2431         }
2432         if (grptp != egrptpPART)
2433         {
2434             if (bVerbose)
2435             {
2436                 fprintf(stderr,
2437                         "Making dummy/rest group for %s containing %d elements\n",
2438                         title, natoms-ntot);
2439             }
2440             /* Add group name "rest" */
2441             grps->nm_ind[grps->nr] = restnm;
2442
2443             /* Assign the rest name to all atoms not currently assigned to a group */
2444             for (j = 0; (j < natoms); j++)
2445             {
2446                 if (cbuf[j] == NOGID)
2447                 {
2448                     cbuf[j] = grps->nr;
2449                 }
2450             }
2451             grps->nr++;
2452         }
2453     }
2454
2455     if (grps->nr == 1 && (ntot == 0 || ntot == natoms))
2456     {
2457         /* All atoms are part of one (or no) group, no index required */
2458         groups->ngrpnr[gtype] = 0;
2459         groups->grpnr[gtype]  = NULL;
2460     }
2461     else
2462     {
2463         groups->ngrpnr[gtype] = natoms;
2464         snew(groups->grpnr[gtype], natoms);
2465         for (j = 0; (j < natoms); j++)
2466         {
2467             groups->grpnr[gtype][j] = cbuf[j];
2468         }
2469     }
2470
2471     sfree(cbuf);
2472
2473     return (bRest && grptp == egrptpPART);
2474 }
2475
2476 static void calc_nrdf(gmx_mtop_t *mtop, t_inputrec *ir, char **gnames)
2477 {
2478     t_grpopts              *opts;
2479     gmx_groups_t           *groups;
2480     t_pull                 *pull;
2481     int                     natoms, ai, aj, i, j, d, g, imin, jmin, nc;
2482     t_iatom                *ia;
2483     int                    *nrdf2, *na_vcm, na_tot;
2484     double                 *nrdf_tc, *nrdf_vcm, nrdf_uc, n_sub = 0;
2485     gmx_mtop_atomloop_all_t aloop;
2486     t_atom                 *atom;
2487     int                     mb, mol, ftype, as;
2488     gmx_molblock_t         *molb;
2489     gmx_moltype_t          *molt;
2490
2491     /* Calculate nrdf.
2492      * First calc 3xnr-atoms for each group
2493      * then subtract half a degree of freedom for each constraint
2494      *
2495      * Only atoms and nuclei contribute to the degrees of freedom...
2496      */
2497
2498     opts = &ir->opts;
2499
2500     groups = &mtop->groups;
2501     natoms = mtop->natoms;
2502
2503     /* Allocate one more for a possible rest group */
2504     /* We need to sum degrees of freedom into doubles,
2505      * since floats give too low nrdf's above 3 million atoms.
2506      */
2507     snew(nrdf_tc, groups->grps[egcTC].nr+1);
2508     snew(nrdf_vcm, groups->grps[egcVCM].nr+1);
2509     snew(na_vcm, groups->grps[egcVCM].nr+1);
2510
2511     for (i = 0; i < groups->grps[egcTC].nr; i++)
2512     {
2513         nrdf_tc[i] = 0;
2514     }
2515     for (i = 0; i < groups->grps[egcVCM].nr+1; i++)
2516     {
2517         nrdf_vcm[i] = 0;
2518     }
2519
2520     snew(nrdf2, natoms);
2521     aloop = gmx_mtop_atomloop_all_init(mtop);
2522     while (gmx_mtop_atomloop_all_next(aloop, &i, &atom))
2523     {
2524         nrdf2[i] = 0;
2525         if (atom->ptype == eptAtom || atom->ptype == eptNucleus)
2526         {
2527             g = ggrpnr(groups, egcFREEZE, i);
2528             /* Double count nrdf for particle i */
2529             for (d = 0; d < DIM; d++)
2530             {
2531                 if (opts->nFreeze[g][d] == 0)
2532                 {
2533                     nrdf2[i] += 2;
2534                 }
2535             }
2536             nrdf_tc [ggrpnr(groups, egcTC, i)]  += 0.5*nrdf2[i];
2537             nrdf_vcm[ggrpnr(groups, egcVCM, i)] += 0.5*nrdf2[i];
2538         }
2539     }
2540
2541     as = 0;
2542     for (mb = 0; mb < mtop->nmolblock; mb++)
2543     {
2544         molb = &mtop->molblock[mb];
2545         molt = &mtop->moltype[molb->type];
2546         atom = molt->atoms.atom;
2547         for (mol = 0; mol < molb->nmol; mol++)
2548         {
2549             for (ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++)
2550             {
2551                 ia = molt->ilist[ftype].iatoms;
2552                 for (i = 0; i < molt->ilist[ftype].nr; )
2553                 {
2554                     /* Subtract degrees of freedom for the constraints,
2555                      * if the particles still have degrees of freedom left.
2556                      * If one of the particles is a vsite or a shell, then all
2557                      * constraint motion will go there, but since they do not
2558                      * contribute to the constraints the degrees of freedom do not
2559                      * change.
2560                      */
2561                     ai = as + ia[1];
2562                     aj = as + ia[2];
2563                     if (((atom[ia[1]].ptype == eptNucleus) ||
2564                          (atom[ia[1]].ptype == eptAtom)) &&
2565                         ((atom[ia[2]].ptype == eptNucleus) ||
2566                          (atom[ia[2]].ptype == eptAtom)))
2567                     {
2568                         if (nrdf2[ai] > 0)
2569                         {
2570                             jmin = 1;
2571                         }
2572                         else
2573                         {
2574                             jmin = 2;
2575                         }
2576                         if (nrdf2[aj] > 0)
2577                         {
2578                             imin = 1;
2579                         }
2580                         else
2581                         {
2582                             imin = 2;
2583                         }
2584                         imin       = min(imin, nrdf2[ai]);
2585                         jmin       = min(jmin, nrdf2[aj]);
2586                         nrdf2[ai] -= imin;
2587                         nrdf2[aj] -= jmin;
2588                         nrdf_tc [ggrpnr(groups, egcTC, ai)]  -= 0.5*imin;
2589                         nrdf_tc [ggrpnr(groups, egcTC, aj)]  -= 0.5*jmin;
2590                         nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5*imin;
2591                         nrdf_vcm[ggrpnr(groups, egcVCM, aj)] -= 0.5*jmin;
2592                     }
2593                     ia += interaction_function[ftype].nratoms+1;
2594                     i  += interaction_function[ftype].nratoms+1;
2595                 }
2596             }
2597             ia = molt->ilist[F_SETTLE].iatoms;
2598             for (i = 0; i < molt->ilist[F_SETTLE].nr; )
2599             {
2600                 /* Subtract 1 dof from every atom in the SETTLE */
2601                 for (j = 0; j < 3; j++)
2602                 {
2603                     ai         = as + ia[1+j];
2604                     imin       = min(2, nrdf2[ai]);
2605                     nrdf2[ai] -= imin;
2606                     nrdf_tc [ggrpnr(groups, egcTC, ai)]  -= 0.5*imin;
2607                     nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5*imin;
2608                 }
2609                 ia += 4;
2610                 i  += 4;
2611             }
2612             as += molt->atoms.nr;
2613         }
2614     }
2615
2616     if (ir->ePull == epullCONSTRAINT)
2617     {
2618         /* Correct nrdf for the COM constraints.
2619          * We correct using the TC and VCM group of the first atom
2620          * in the reference and pull group. If atoms in one pull group
2621          * belong to different TC or VCM groups it is anyhow difficult
2622          * to determine the optimal nrdf assignment.
2623          */
2624         pull = ir->pull;
2625         if (pull->eGeom == epullgPOS)
2626         {
2627             nc = 0;
2628             for (i = 0; i < DIM; i++)
2629             {
2630                 if (pull->dim[i])
2631                 {
2632                     nc++;
2633                 }
2634             }
2635         }
2636         else
2637         {
2638             nc = 1;
2639         }
2640         for (i = 0; i < pull->ngrp; i++)
2641         {
2642             imin = 2*nc;
2643             if (pull->grp[0].nat > 0)
2644             {
2645                 /* Subtract 1/2 dof from the reference group */
2646                 ai = pull->grp[0].ind[0];
2647                 if (nrdf_tc[ggrpnr(groups, egcTC, ai)] > 1)
2648                 {
2649                     nrdf_tc [ggrpnr(groups, egcTC, ai)]  -= 0.5;
2650                     nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5;
2651                     imin--;
2652                 }
2653             }
2654             /* Subtract 1/2 dof from the pulled group */
2655             ai = pull->grp[1+i].ind[0];
2656             nrdf_tc [ggrpnr(groups, egcTC, ai)]  -= 0.5*imin;
2657             nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5*imin;
2658             if (nrdf_tc[ggrpnr(groups, egcTC, ai)] < 0)
2659             {
2660                 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)]]);
2661             }
2662         }
2663     }
2664
2665     if (ir->nstcomm != 0)
2666     {
2667         /* Subtract 3 from the number of degrees of freedom in each vcm group
2668          * when com translation is removed and 6 when rotation is removed
2669          * as well.
2670          */
2671         switch (ir->comm_mode)
2672         {
2673             case ecmLINEAR:
2674                 n_sub = ndof_com(ir);
2675                 break;
2676             case ecmANGULAR:
2677                 n_sub = 6;
2678                 break;
2679             default:
2680                 n_sub = 0;
2681                 gmx_incons("Checking comm_mode");
2682         }
2683
2684         for (i = 0; i < groups->grps[egcTC].nr; i++)
2685         {
2686             /* Count the number of atoms of TC group i for every VCM group */
2687             for (j = 0; j < groups->grps[egcVCM].nr+1; j++)
2688             {
2689                 na_vcm[j] = 0;
2690             }
2691             na_tot = 0;
2692             for (ai = 0; ai < natoms; ai++)
2693             {
2694                 if (ggrpnr(groups, egcTC, ai) == i)
2695                 {
2696                     na_vcm[ggrpnr(groups, egcVCM, ai)]++;
2697                     na_tot++;
2698                 }
2699             }
2700             /* Correct for VCM removal according to the fraction of each VCM
2701              * group present in this TC group.
2702              */
2703             nrdf_uc = nrdf_tc[i];
2704             if (debug)
2705             {
2706                 fprintf(debug, "T-group[%d] nrdf_uc = %g, n_sub = %g\n",
2707                         i, nrdf_uc, n_sub);
2708             }
2709             nrdf_tc[i] = 0;
2710             for (j = 0; j < groups->grps[egcVCM].nr+1; j++)
2711             {
2712                 if (nrdf_vcm[j] > n_sub)
2713                 {
2714                     nrdf_tc[i] += nrdf_uc*((double)na_vcm[j]/(double)na_tot)*
2715                         (nrdf_vcm[j] - n_sub)/nrdf_vcm[j];
2716                 }
2717                 if (debug)
2718                 {
2719                     fprintf(debug, "  nrdf_vcm[%d] = %g, nrdf = %g\n",
2720                             j, nrdf_vcm[j], nrdf_tc[i]);
2721                 }
2722             }
2723         }
2724     }
2725     for (i = 0; (i < groups->grps[egcTC].nr); i++)
2726     {
2727         opts->nrdf[i] = nrdf_tc[i];
2728         if (opts->nrdf[i] < 0)
2729         {
2730             opts->nrdf[i] = 0;
2731         }
2732         fprintf(stderr,
2733                 "Number of degrees of freedom in T-Coupling group %s is %.2f\n",
2734                 gnames[groups->grps[egcTC].nm_ind[i]], opts->nrdf[i]);
2735     }
2736
2737     sfree(nrdf2);
2738     sfree(nrdf_tc);
2739     sfree(nrdf_vcm);
2740     sfree(na_vcm);
2741 }
2742
2743 static void decode_cos(char *s, t_cosines *cosine)
2744 {
2745     char   *t;
2746     char    format[STRLEN], f1[STRLEN];
2747     double  a, phi;
2748     int     i;
2749
2750     t = strdup(s);
2751     trim(t);
2752
2753     cosine->n   = 0;
2754     cosine->a   = NULL;
2755     cosine->phi = NULL;
2756     if (strlen(t))
2757     {
2758         sscanf(t, "%d", &(cosine->n));
2759         if (cosine->n <= 0)
2760         {
2761             cosine->n = 0;
2762         }
2763         else
2764         {
2765             snew(cosine->a, cosine->n);
2766             snew(cosine->phi, cosine->n);
2767
2768             sprintf(format, "%%*d");
2769             for (i = 0; (i < cosine->n); i++)
2770             {
2771                 strcpy(f1, format);
2772                 strcat(f1, "%lf%lf");
2773                 if (sscanf(t, f1, &a, &phi) < 2)
2774                 {
2775                     gmx_fatal(FARGS, "Invalid input for electric field shift: '%s'", t);
2776                 }
2777                 cosine->a[i]   = a;
2778                 cosine->phi[i] = phi;
2779                 strcat(format, "%*lf%*lf");
2780             }
2781         }
2782     }
2783     sfree(t);
2784 }
2785
2786 static gmx_bool do_egp_flag(t_inputrec *ir, gmx_groups_t *groups,
2787                             const char *option, const char *val, int flag)
2788 {
2789     /* The maximum number of energy group pairs would be MAXPTR*(MAXPTR+1)/2.
2790      * But since this is much larger than STRLEN, such a line can not be parsed.
2791      * The real maximum is the number of names that fit in a string: STRLEN/2.
2792      */
2793 #define EGP_MAX (STRLEN/2)
2794     int      nelem, i, j, k, nr;
2795     char    *names[EGP_MAX];
2796     char  ***gnames;
2797     gmx_bool bSet;
2798
2799     gnames = groups->grpname;
2800
2801     nelem = str_nelem(val, EGP_MAX, names);
2802     if (nelem % 2 != 0)
2803     {
2804         gmx_fatal(FARGS, "The number of groups for %s is odd", option);
2805     }
2806     nr   = groups->grps[egcENER].nr;
2807     bSet = FALSE;
2808     for (i = 0; i < nelem/2; i++)
2809     {
2810         j = 0;
2811         while ((j < nr) &&
2812                gmx_strcasecmp(names[2*i], *(gnames[groups->grps[egcENER].nm_ind[j]])))
2813         {
2814             j++;
2815         }
2816         if (j == nr)
2817         {
2818             gmx_fatal(FARGS, "%s in %s is not an energy group\n",
2819                       names[2*i], option);
2820         }
2821         k = 0;
2822         while ((k < nr) &&
2823                gmx_strcasecmp(names[2*i+1], *(gnames[groups->grps[egcENER].nm_ind[k]])))
2824         {
2825             k++;
2826         }
2827         if (k == nr)
2828         {
2829             gmx_fatal(FARGS, "%s in %s is not an energy group\n",
2830                       names[2*i+1], option);
2831         }
2832         if ((j < nr) && (k < nr))
2833         {
2834             ir->opts.egp_flags[nr*j+k] |= flag;
2835             ir->opts.egp_flags[nr*k+j] |= flag;
2836             bSet = TRUE;
2837         }
2838     }
2839
2840     return bSet;
2841 }
2842
2843 void do_index(const char* mdparin, const char *ndx,
2844               gmx_mtop_t *mtop,
2845               gmx_bool bVerbose,
2846               t_inputrec *ir, rvec *v,
2847               warninp_t wi)
2848 {
2849     t_blocka     *grps;
2850     gmx_groups_t *groups;
2851     int           natoms;
2852     t_symtab     *symtab;
2853     t_atoms       atoms_all;
2854     char          warnbuf[STRLEN], **gnames;
2855     int           nr, ntcg, ntau_t, nref_t, nacc, nofg, nSA, nSA_points, nSA_time, nSA_temp;
2856     real          tau_min;
2857     int           nstcmin;
2858     int           nacg, nfreeze, nfrdim, nenergy, nvcm, nuser;
2859     char         *ptr1[MAXPTR], *ptr2[MAXPTR], *ptr3[MAXPTR];
2860     int           i, j, k, restnm;
2861     real          SAtime;
2862     gmx_bool      bExcl, bTable, bSetTCpar, bAnneal, bRest;
2863     int           nQMmethod, nQMbasis, nQMcharge, nQMmult, nbSH, nCASorb, nCASelec,
2864                   nSAon, nSAoff, nSAsteps, nQMg, nbOPT, nbTS;
2865     char          warn_buf[STRLEN];
2866
2867     if (bVerbose)
2868     {
2869         fprintf(stderr, "processing index file...\n");
2870     }
2871     debug_gmx();
2872     if (ndx == NULL)
2873     {
2874         snew(grps, 1);
2875         snew(grps->index, 1);
2876         snew(gnames, 1);
2877         atoms_all = gmx_mtop_global_atoms(mtop);
2878         analyse(&atoms_all, grps, &gnames, FALSE, TRUE);
2879         free_t_atoms(&atoms_all, FALSE);
2880     }
2881     else
2882     {
2883         grps = init_index(ndx, &gnames);
2884     }
2885
2886     groups = &mtop->groups;
2887     natoms = mtop->natoms;
2888     symtab = &mtop->symtab;
2889
2890     snew(groups->grpname, grps->nr+1);
2891
2892     for (i = 0; (i < grps->nr); i++)
2893     {
2894         groups->grpname[i] = put_symtab(symtab, gnames[i]);
2895     }
2896     groups->grpname[i] = put_symtab(symtab, "rest");
2897     restnm             = i;
2898     srenew(gnames, grps->nr+1);
2899     gnames[restnm]   = *(groups->grpname[i]);
2900     groups->ngrpname = grps->nr+1;
2901
2902     set_warning_line(wi, mdparin, -1);
2903
2904     ntau_t = str_nelem(tau_t, MAXPTR, ptr1);
2905     nref_t = str_nelem(ref_t, MAXPTR, ptr2);
2906     ntcg   = str_nelem(tcgrps, MAXPTR, ptr3);
2907     if ((ntau_t != ntcg) || (nref_t != ntcg))
2908     {
2909         gmx_fatal(FARGS, "Invalid T coupling input: %d groups, %d ref-t values and "
2910                   "%d tau-t values", ntcg, nref_t, ntau_t);
2911     }
2912
2913     bSetTCpar = (ir->etc || EI_SD(ir->eI) || ir->eI == eiBD || EI_TPI(ir->eI));
2914     do_numbering(natoms, groups, ntcg, ptr3, grps, gnames, egcTC,
2915                  restnm, bSetTCpar ? egrptpALL : egrptpALL_GENREST, bVerbose, wi);
2916     nr            = groups->grps[egcTC].nr;
2917     ir->opts.ngtc = nr;
2918     snew(ir->opts.nrdf, nr);
2919     snew(ir->opts.tau_t, nr);
2920     snew(ir->opts.ref_t, nr);
2921     if (ir->eI == eiBD && ir->bd_fric == 0)
2922     {
2923         fprintf(stderr, "bd-fric=0, so tau-t will be used as the inverse friction constant(s)\n");
2924     }
2925
2926     if (bSetTCpar)
2927     {
2928         if (nr != nref_t)
2929         {
2930             gmx_fatal(FARGS, "Not enough ref-t and tau-t values!");
2931         }
2932
2933         tau_min = 1e20;
2934         for (i = 0; (i < nr); i++)
2935         {
2936             ir->opts.tau_t[i] = strtod(ptr1[i], NULL);
2937             if ((ir->eI == eiBD || ir->eI == eiSD2) && ir->opts.tau_t[i] <= 0)
2938             {
2939                 sprintf(warn_buf, "With integrator %s tau-t should be larger than 0", ei_names[ir->eI]);
2940                 warning_error(wi, warn_buf);
2941             }
2942
2943             if (ir->etc != etcVRESCALE && ir->opts.tau_t[i] == 0)
2944             {
2945                 warning_note(wi, "tau-t = -1 is the value to signal that a group should not have temperature coupling. Treating your use of tau-t = 0 as if you used -1.");
2946             }
2947
2948             if (ir->opts.tau_t[i] >= 0)
2949             {
2950                 tau_min = min(tau_min, ir->opts.tau_t[i]);
2951             }
2952         }
2953         if (ir->etc != etcNO && ir->nsttcouple == -1)
2954         {
2955             ir->nsttcouple = ir_optimal_nsttcouple(ir);
2956         }
2957
2958         if (EI_VV(ir->eI))
2959         {
2960             if ((ir->etc == etcNOSEHOOVER) && (ir->epc == epcBERENDSEN))
2961             {
2962                 gmx_fatal(FARGS, "Cannot do Nose-Hoover temperature with Berendsen pressure control with md-vv; use either vrescale temperature with berendsen pressure or Nose-Hoover temperature with MTTK pressure");
2963             }
2964             if ((ir->epc == epcMTTK) && (ir->etc > etcNO))
2965             {
2966                 if (ir->nstpcouple != ir->nsttcouple)
2967                 {
2968                     int mincouple = min(ir->nstpcouple, ir->nsttcouple);
2969                     ir->nstpcouple = ir->nsttcouple = mincouple;
2970                     sprintf(warn_buf, "for current Trotter decomposition methods with vv, nsttcouple and nstpcouple must be equal.  Both have been reset to min(nsttcouple,nstpcouple) = %d", mincouple);
2971                     warning_note(wi, warn_buf);
2972                 }
2973             }
2974         }
2975         /* velocity verlet with averaged kinetic energy KE = 0.5*(v(t+1/2) - v(t-1/2)) is implemented
2976            primarily for testing purposes, and does not work with temperature coupling other than 1 */
2977
2978         if (ETC_ANDERSEN(ir->etc))
2979         {
2980             if (ir->nsttcouple != 1)
2981             {
2982                 ir->nsttcouple = 1;
2983                 sprintf(warn_buf, "Andersen temperature control methods assume nsttcouple = 1; there is no need for larger nsttcouple > 1, since no global parameters are computed. nsttcouple has been reset to 1");
2984                 warning_note(wi, warn_buf);
2985             }
2986         }
2987         nstcmin = tcouple_min_integration_steps(ir->etc);
2988         if (nstcmin > 1)
2989         {
2990             if (tau_min/(ir->delta_t*ir->nsttcouple) < nstcmin)
2991             {
2992                 sprintf(warn_buf, "For proper integration of the %s thermostat, tau-t (%g) should be at least %d times larger than nsttcouple*dt (%g)",
2993                         ETCOUPLTYPE(ir->etc),
2994                         tau_min, nstcmin,
2995                         ir->nsttcouple*ir->delta_t);
2996                 warning(wi, warn_buf);
2997             }
2998         }
2999         for (i = 0; (i < nr); i++)
3000         {
3001             ir->opts.ref_t[i] = strtod(ptr2[i], NULL);
3002             if (ir->opts.ref_t[i] < 0)
3003             {
3004                 gmx_fatal(FARGS, "ref-t for group %d negative", i);
3005             }
3006         }
3007         /* set the lambda mc temperature to the md integrator temperature (which should be defined
3008            if we are in this conditional) if mc_temp is negative */
3009         if (ir->expandedvals->mc_temp < 0)
3010         {
3011             ir->expandedvals->mc_temp = ir->opts.ref_t[0]; /*for now, set to the first reft */
3012         }
3013     }
3014
3015     /* Simulated annealing for each group. There are nr groups */
3016     nSA = str_nelem(anneal, MAXPTR, ptr1);
3017     if (nSA == 1 && (ptr1[0][0] == 'n' || ptr1[0][0] == 'N'))
3018     {
3019         nSA = 0;
3020     }
3021     if (nSA > 0 && nSA != nr)
3022     {
3023         gmx_fatal(FARGS, "Not enough annealing values: %d (for %d groups)\n", nSA, nr);
3024     }
3025     else
3026     {
3027         snew(ir->opts.annealing, nr);
3028         snew(ir->opts.anneal_npoints, nr);
3029         snew(ir->opts.anneal_time, nr);
3030         snew(ir->opts.anneal_temp, nr);
3031         for (i = 0; i < nr; i++)
3032         {
3033             ir->opts.annealing[i]      = eannNO;
3034             ir->opts.anneal_npoints[i] = 0;
3035             ir->opts.anneal_time[i]    = NULL;
3036             ir->opts.anneal_temp[i]    = NULL;
3037         }
3038         if (nSA > 0)
3039         {
3040             bAnneal = FALSE;
3041             for (i = 0; i < nr; i++)
3042             {
3043                 if (ptr1[i][0] == 'n' || ptr1[i][0] == 'N')
3044                 {
3045                     ir->opts.annealing[i] = eannNO;
3046                 }
3047                 else if (ptr1[i][0] == 's' || ptr1[i][0] == 'S')
3048                 {
3049                     ir->opts.annealing[i] = eannSINGLE;
3050                     bAnneal               = TRUE;
3051                 }
3052                 else if (ptr1[i][0] == 'p' || ptr1[i][0] == 'P')
3053                 {
3054                     ir->opts.annealing[i] = eannPERIODIC;
3055                     bAnneal               = TRUE;
3056                 }
3057             }
3058             if (bAnneal)
3059             {
3060                 /* Read the other fields too */
3061                 nSA_points = str_nelem(anneal_npoints, MAXPTR, ptr1);
3062                 if (nSA_points != nSA)
3063                 {
3064                     gmx_fatal(FARGS, "Found %d annealing-npoints values for %d groups\n", nSA_points, nSA);
3065                 }
3066                 for (k = 0, i = 0; i < nr; i++)
3067                 {
3068                     ir->opts.anneal_npoints[i] = strtol(ptr1[i], NULL, 10);
3069                     if (ir->opts.anneal_npoints[i] == 1)
3070                     {
3071                         gmx_fatal(FARGS, "Please specify at least a start and an end point for annealing\n");
3072                     }
3073                     snew(ir->opts.anneal_time[i], ir->opts.anneal_npoints[i]);
3074                     snew(ir->opts.anneal_temp[i], ir->opts.anneal_npoints[i]);
3075                     k += ir->opts.anneal_npoints[i];
3076                 }
3077
3078                 nSA_time = str_nelem(anneal_time, MAXPTR, ptr1);
3079                 if (nSA_time != k)
3080                 {
3081                     gmx_fatal(FARGS, "Found %d annealing-time values, wanter %d\n", nSA_time, k);
3082                 }
3083                 nSA_temp = str_nelem(anneal_temp, MAXPTR, ptr2);
3084                 if (nSA_temp != k)
3085                 {
3086                     gmx_fatal(FARGS, "Found %d annealing-temp values, wanted %d\n", nSA_temp, k);
3087                 }
3088
3089                 for (i = 0, k = 0; i < nr; i++)
3090                 {
3091
3092                     for (j = 0; j < ir->opts.anneal_npoints[i]; j++)
3093                     {
3094                         ir->opts.anneal_time[i][j] = strtod(ptr1[k], NULL);
3095                         ir->opts.anneal_temp[i][j] = strtod(ptr2[k], NULL);
3096                         if (j == 0)
3097                         {
3098                             if (ir->opts.anneal_time[i][0] > (ir->init_t+GMX_REAL_EPS))
3099                             {
3100                                 gmx_fatal(FARGS, "First time point for annealing > init_t.\n");
3101                             }
3102                         }
3103                         else
3104                         {
3105                             /* j>0 */
3106                             if (ir->opts.anneal_time[i][j] < ir->opts.anneal_time[i][j-1])
3107                             {
3108                                 gmx_fatal(FARGS, "Annealing timepoints out of order: t=%f comes after t=%f\n",
3109                                           ir->opts.anneal_time[i][j], ir->opts.anneal_time[i][j-1]);
3110                             }
3111                         }
3112                         if (ir->opts.anneal_temp[i][j] < 0)
3113                         {
3114                             gmx_fatal(FARGS, "Found negative temperature in annealing: %f\n", ir->opts.anneal_temp[i][j]);
3115                         }
3116                         k++;
3117                     }
3118                 }
3119                 /* Print out some summary information, to make sure we got it right */
3120                 for (i = 0, k = 0; i < nr; i++)
3121                 {
3122                     if (ir->opts.annealing[i] != eannNO)
3123                     {
3124                         j = groups->grps[egcTC].nm_ind[i];
3125                         fprintf(stderr, "Simulated annealing for group %s: %s, %d timepoints\n",
3126                                 *(groups->grpname[j]), eann_names[ir->opts.annealing[i]],
3127                                 ir->opts.anneal_npoints[i]);
3128                         fprintf(stderr, "Time (ps)   Temperature (K)\n");
3129                         /* All terms except the last one */
3130                         for (j = 0; j < (ir->opts.anneal_npoints[i]-1); j++)
3131                         {
3132                             fprintf(stderr, "%9.1f      %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3133                         }
3134
3135                         /* Finally the last one */
3136                         j = ir->opts.anneal_npoints[i]-1;
3137                         if (ir->opts.annealing[i] == eannSINGLE)
3138                         {
3139                             fprintf(stderr, "%9.1f-     %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3140                         }
3141                         else
3142                         {
3143                             fprintf(stderr, "%9.1f      %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3144                             if (fabs(ir->opts.anneal_temp[i][j]-ir->opts.anneal_temp[i][0]) > GMX_REAL_EPS)
3145                             {
3146                                 warning_note(wi, "There is a temperature jump when your annealing loops back.\n");
3147                             }
3148                         }
3149                     }
3150                 }
3151             }
3152         }
3153     }
3154
3155     if (ir->ePull != epullNO)
3156     {
3157         make_pull_groups(ir->pull, pull_grp, grps, gnames);
3158     }
3159
3160     if (ir->bRot)
3161     {
3162         make_rotation_groups(ir->rot, rot_grp, grps, gnames);
3163     }
3164
3165     nacc = str_nelem(acc, MAXPTR, ptr1);
3166     nacg = str_nelem(accgrps, MAXPTR, ptr2);
3167     if (nacg*DIM != nacc)
3168     {
3169         gmx_fatal(FARGS, "Invalid Acceleration input: %d groups and %d acc. values",
3170                   nacg, nacc);
3171     }
3172     do_numbering(natoms, groups, nacg, ptr2, grps, gnames, egcACC,
3173                  restnm, egrptpALL_GENREST, bVerbose, wi);
3174     nr = groups->grps[egcACC].nr;
3175     snew(ir->opts.acc, nr);
3176     ir->opts.ngacc = nr;
3177
3178     for (i = k = 0; (i < nacg); i++)
3179     {
3180         for (j = 0; (j < DIM); j++, k++)
3181         {
3182             ir->opts.acc[i][j] = strtod(ptr1[k], NULL);
3183         }
3184     }
3185     for (; (i < nr); i++)
3186     {
3187         for (j = 0; (j < DIM); j++)
3188         {
3189             ir->opts.acc[i][j] = 0;
3190         }
3191     }
3192
3193     nfrdim  = str_nelem(frdim, MAXPTR, ptr1);
3194     nfreeze = str_nelem(freeze, MAXPTR, ptr2);
3195     if (nfrdim != DIM*nfreeze)
3196     {
3197         gmx_fatal(FARGS, "Invalid Freezing input: %d groups and %d freeze values",
3198                   nfreeze, nfrdim);
3199     }
3200     do_numbering(natoms, groups, nfreeze, ptr2, grps, gnames, egcFREEZE,
3201                  restnm, egrptpALL_GENREST, bVerbose, wi);
3202     nr             = groups->grps[egcFREEZE].nr;
3203     ir->opts.ngfrz = nr;
3204     snew(ir->opts.nFreeze, nr);
3205     for (i = k = 0; (i < nfreeze); i++)
3206     {
3207         for (j = 0; (j < DIM); j++, k++)
3208         {
3209             ir->opts.nFreeze[i][j] = (gmx_strncasecmp(ptr1[k], "Y", 1) == 0);
3210             if (!ir->opts.nFreeze[i][j])
3211             {
3212                 if (gmx_strncasecmp(ptr1[k], "N", 1) != 0)
3213                 {
3214                     sprintf(warnbuf, "Please use Y(ES) or N(O) for freezedim only "
3215                             "(not %s)", ptr1[k]);
3216                     warning(wi, warn_buf);
3217                 }
3218             }
3219         }
3220     }
3221     for (; (i < nr); i++)
3222     {
3223         for (j = 0; (j < DIM); j++)
3224         {
3225             ir->opts.nFreeze[i][j] = 0;
3226         }
3227     }
3228
3229     nenergy = str_nelem(energy, MAXPTR, ptr1);
3230     do_numbering(natoms, groups, nenergy, ptr1, grps, gnames, egcENER,
3231                  restnm, egrptpALL_GENREST, bVerbose, wi);
3232     add_wall_energrps(groups, ir->nwall, symtab);
3233     ir->opts.ngener = groups->grps[egcENER].nr;
3234     nvcm            = str_nelem(vcm, MAXPTR, ptr1);
3235     bRest           =
3236         do_numbering(natoms, groups, nvcm, ptr1, grps, gnames, egcVCM,
3237                      restnm, nvcm == 0 ? egrptpALL_GENREST : egrptpPART, bVerbose, wi);
3238     if (bRest)
3239     {
3240         warning(wi, "Some atoms are not part of any center of mass motion removal group.\n"
3241                 "This may lead to artifacts.\n"
3242                 "In most cases one should use one group for the whole system.");
3243     }
3244
3245     /* Now we have filled the freeze struct, so we can calculate NRDF */
3246     calc_nrdf(mtop, ir, gnames);
3247
3248     if (v && NULL)
3249     {
3250         real fac, ntot = 0;
3251
3252         /* Must check per group! */
3253         for (i = 0; (i < ir->opts.ngtc); i++)
3254         {
3255             ntot += ir->opts.nrdf[i];
3256         }
3257         if (ntot != (DIM*natoms))
3258         {
3259             fac = sqrt(ntot/(DIM*natoms));
3260             if (bVerbose)
3261             {
3262                 fprintf(stderr, "Scaling velocities by a factor of %.3f to account for constraints\n"
3263                         "and removal of center of mass motion\n", fac);
3264             }
3265             for (i = 0; (i < natoms); i++)
3266             {
3267                 svmul(fac, v[i], v[i]);
3268             }
3269         }
3270     }
3271
3272     nuser = str_nelem(user1, MAXPTR, ptr1);
3273     do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcUser1,
3274                  restnm, egrptpALL_GENREST, bVerbose, wi);
3275     nuser = str_nelem(user2, MAXPTR, ptr1);
3276     do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcUser2,
3277                  restnm, egrptpALL_GENREST, bVerbose, wi);
3278     nuser = str_nelem(xtc_grps, MAXPTR, ptr1);
3279     do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcXTC,
3280                  restnm, egrptpONE, bVerbose, wi);
3281     nofg = str_nelem(orirefitgrp, MAXPTR, ptr1);
3282     do_numbering(natoms, groups, nofg, ptr1, grps, gnames, egcORFIT,
3283                  restnm, egrptpALL_GENREST, bVerbose, wi);
3284
3285     /* QMMM input processing */
3286     nQMg          = str_nelem(QMMM, MAXPTR, ptr1);
3287     nQMmethod     = str_nelem(QMmethod, MAXPTR, ptr2);
3288     nQMbasis      = str_nelem(QMbasis, MAXPTR, ptr3);
3289     if ((nQMmethod != nQMg) || (nQMbasis != nQMg))
3290     {
3291         gmx_fatal(FARGS, "Invalid QMMM input: %d groups %d basissets"
3292                   " and %d methods\n", nQMg, nQMbasis, nQMmethod);
3293     }
3294     /* group rest, if any, is always MM! */
3295     do_numbering(natoms, groups, nQMg, ptr1, grps, gnames, egcQMMM,
3296                  restnm, egrptpALL_GENREST, bVerbose, wi);
3297     nr            = nQMg; /*atoms->grps[egcQMMM].nr;*/
3298     ir->opts.ngQM = nQMg;
3299     snew(ir->opts.QMmethod, nr);
3300     snew(ir->opts.QMbasis, nr);
3301     for (i = 0; i < nr; i++)
3302     {
3303         /* input consists of strings: RHF CASSCF PM3 .. These need to be
3304          * converted to the corresponding enum in names.c
3305          */
3306         ir->opts.QMmethod[i] = search_QMstring(ptr2[i], eQMmethodNR,
3307                                                eQMmethod_names);
3308         ir->opts.QMbasis[i]  = search_QMstring(ptr3[i], eQMbasisNR,
3309                                                eQMbasis_names);
3310
3311     }
3312     nQMmult   = str_nelem(QMmult, MAXPTR, ptr1);
3313     nQMcharge = str_nelem(QMcharge, MAXPTR, ptr2);
3314     nbSH      = str_nelem(bSH, MAXPTR, ptr3);
3315     snew(ir->opts.QMmult, nr);
3316     snew(ir->opts.QMcharge, nr);
3317     snew(ir->opts.bSH, nr);
3318
3319     for (i = 0; i < nr; i++)
3320     {
3321         ir->opts.QMmult[i]   = strtol(ptr1[i], NULL, 10);
3322         ir->opts.QMcharge[i] = strtol(ptr2[i], NULL, 10);
3323         ir->opts.bSH[i]      = (gmx_strncasecmp(ptr3[i], "Y", 1) == 0);
3324     }
3325
3326     nCASelec  = str_nelem(CASelectrons, MAXPTR, ptr1);
3327     nCASorb   = str_nelem(CASorbitals, MAXPTR, ptr2);
3328     snew(ir->opts.CASelectrons, nr);
3329     snew(ir->opts.CASorbitals, nr);
3330     for (i = 0; i < nr; i++)
3331     {
3332         ir->opts.CASelectrons[i] = strtol(ptr1[i], NULL, 10);
3333         ir->opts.CASorbitals[i]  = strtol(ptr2[i], NULL, 10);
3334     }
3335     /* special optimization options */
3336
3337     nbOPT = str_nelem(bOPT, MAXPTR, ptr1);
3338     nbTS  = str_nelem(bTS, MAXPTR, ptr2);
3339     snew(ir->opts.bOPT, nr);
3340     snew(ir->opts.bTS, nr);
3341     for (i = 0; i < nr; i++)
3342     {
3343         ir->opts.bOPT[i] = (gmx_strncasecmp(ptr1[i], "Y", 1) == 0);
3344         ir->opts.bTS[i]  = (gmx_strncasecmp(ptr2[i], "Y", 1) == 0);
3345     }
3346     nSAon     = str_nelem(SAon, MAXPTR, ptr1);
3347     nSAoff    = str_nelem(SAoff, MAXPTR, ptr2);
3348     nSAsteps  = str_nelem(SAsteps, MAXPTR, ptr3);
3349     snew(ir->opts.SAon, nr);
3350     snew(ir->opts.SAoff, nr);
3351     snew(ir->opts.SAsteps, nr);
3352
3353     for (i = 0; i < nr; i++)
3354     {
3355         ir->opts.SAon[i]    = strtod(ptr1[i], NULL);
3356         ir->opts.SAoff[i]   = strtod(ptr2[i], NULL);
3357         ir->opts.SAsteps[i] = strtol(ptr3[i], NULL, 10);
3358     }
3359     /* end of QMMM input */
3360
3361     if (bVerbose)
3362     {
3363         for (i = 0; (i < egcNR); i++)
3364         {
3365             fprintf(stderr, "%-16s has %d element(s):", gtypes[i], groups->grps[i].nr);
3366             for (j = 0; (j < groups->grps[i].nr); j++)
3367             {
3368                 fprintf(stderr, " %s", *(groups->grpname[groups->grps[i].nm_ind[j]]));
3369             }
3370             fprintf(stderr, "\n");
3371         }
3372     }
3373
3374     nr = groups->grps[egcENER].nr;
3375     snew(ir->opts.egp_flags, nr*nr);
3376
3377     bExcl = do_egp_flag(ir, groups, "energygrp-excl", egpexcl, EGP_EXCL);
3378     if (bExcl && ir->cutoff_scheme == ecutsVERLET)
3379     {
3380         warning_error(wi, "Energy group exclusions are not (yet) implemented for the Verlet scheme");
3381     }
3382     if (bExcl && EEL_FULL(ir->coulombtype))
3383     {
3384         warning(wi, "Can not exclude the lattice Coulomb energy between energy groups");
3385     }
3386
3387     bTable = do_egp_flag(ir, groups, "energygrp-table", egptable, EGP_TABLE);
3388     if (bTable && !(ir->vdwtype == evdwUSER) &&
3389         !(ir->coulombtype == eelUSER) && !(ir->coulombtype == eelPMEUSER) &&
3390         !(ir->coulombtype == eelPMEUSERSWITCH))
3391     {
3392         gmx_fatal(FARGS, "Can only have energy group pair tables in combination with user tables for VdW and/or Coulomb");
3393     }
3394
3395     decode_cos(efield_x, &(ir->ex[XX]));
3396     decode_cos(efield_xt, &(ir->et[XX]));
3397     decode_cos(efield_y, &(ir->ex[YY]));
3398     decode_cos(efield_yt, &(ir->et[YY]));
3399     decode_cos(efield_z, &(ir->ex[ZZ]));
3400     decode_cos(efield_zt, &(ir->et[ZZ]));
3401
3402     if (ir->bAdress)
3403     {
3404         do_adress_index(ir->adress, groups, gnames, &(ir->opts), wi);
3405     }
3406
3407     for (i = 0; (i < grps->nr); i++)
3408     {
3409         sfree(gnames[i]);
3410     }
3411     sfree(gnames);
3412     done_blocka(grps);
3413     sfree(grps);
3414
3415 }
3416
3417
3418
3419 static void check_disre(gmx_mtop_t *mtop)
3420 {
3421     gmx_ffparams_t *ffparams;
3422     t_functype     *functype;
3423     t_iparams      *ip;
3424     int             i, ndouble, ftype;
3425     int             label, old_label;
3426
3427     if (gmx_mtop_ftype_count(mtop, F_DISRES) > 0)
3428     {
3429         ffparams  = &mtop->ffparams;
3430         functype  = ffparams->functype;
3431         ip        = ffparams->iparams;
3432         ndouble   = 0;
3433         old_label = -1;
3434         for (i = 0; i < ffparams->ntypes; i++)
3435         {
3436             ftype = functype[i];
3437             if (ftype == F_DISRES)
3438             {
3439                 label = ip[i].disres.label;
3440                 if (label == old_label)
3441                 {
3442                     fprintf(stderr, "Distance restraint index %d occurs twice\n", label);
3443                     ndouble++;
3444                 }
3445                 old_label = label;
3446             }
3447         }
3448         if (ndouble > 0)
3449         {
3450             gmx_fatal(FARGS, "Found %d double distance restraint indices,\n"
3451                       "probably the parameters for multiple pairs in one restraint "
3452                       "are not identical\n", ndouble);
3453         }
3454     }
3455 }
3456
3457 static gmx_bool absolute_reference(t_inputrec *ir, gmx_mtop_t *sys,
3458                                    gmx_bool posres_only,
3459                                    ivec AbsRef)
3460 {
3461     int                  d, g, i;
3462     gmx_mtop_ilistloop_t iloop;
3463     t_ilist             *ilist;
3464     int                  nmol;
3465     t_iparams           *pr;
3466
3467     clear_ivec(AbsRef);
3468
3469     if (!posres_only)
3470     {
3471         /* Check the COM */
3472         for (d = 0; d < DIM; d++)
3473         {
3474             AbsRef[d] = (d < ndof_com(ir) ? 0 : 1);
3475         }
3476         /* Check for freeze groups */
3477         for (g = 0; g < ir->opts.ngfrz; g++)
3478         {
3479             for (d = 0; d < DIM; d++)
3480             {
3481                 if (ir->opts.nFreeze[g][d] != 0)
3482                 {
3483                     AbsRef[d] = 1;
3484                 }
3485             }
3486         }
3487     }
3488
3489     /* Check for position restraints */
3490     iloop = gmx_mtop_ilistloop_init(sys);
3491     while (gmx_mtop_ilistloop_next(iloop, &ilist, &nmol))
3492     {
3493         if (nmol > 0 &&
3494             (AbsRef[XX] == 0 || AbsRef[YY] == 0 || AbsRef[ZZ] == 0))
3495         {
3496             for (i = 0; i < ilist[F_POSRES].nr; i += 2)
3497             {
3498                 pr = &sys->ffparams.iparams[ilist[F_POSRES].iatoms[i]];
3499                 for (d = 0; d < DIM; d++)
3500                 {
3501                     if (pr->posres.fcA[d] != 0)
3502                     {
3503                         AbsRef[d] = 1;
3504                     }
3505                 }
3506             }
3507             for (i = 0; i < ilist[F_FBPOSRES].nr; i += 2)
3508             {
3509                 /* Check for flat-bottom posres */
3510                 pr = &sys->ffparams.iparams[ilist[F_FBPOSRES].iatoms[i]];
3511                 if (pr->fbposres.k != 0)
3512                 {
3513                     switch (pr->fbposres.geom)
3514                     {
3515                         case efbposresSPHERE:
3516                             AbsRef[XX] = AbsRef[YY] = AbsRef[ZZ] = 1;
3517                             break;
3518                         case efbposresCYLINDER:
3519                             AbsRef[XX] = AbsRef[YY] = 1;
3520                             break;
3521                         case efbposresX: /* d=XX */
3522                         case efbposresY: /* d=YY */
3523                         case efbposresZ: /* d=ZZ */
3524                             d         = pr->fbposres.geom - efbposresX;
3525                             AbsRef[d] = 1;
3526                             break;
3527                         default:
3528                             gmx_fatal(FARGS, " Invalid geometry for flat-bottom position restraint.\n"
3529                                       "Expected nr between 1 and %d. Found %d\n", efbposresNR-1,
3530                                       pr->fbposres.geom);
3531                     }
3532                 }
3533             }
3534         }
3535     }
3536
3537     return (AbsRef[XX] != 0 && AbsRef[YY] != 0 && AbsRef[ZZ] != 0);
3538 }
3539
3540 void triple_check(const char *mdparin, t_inputrec *ir, gmx_mtop_t *sys,
3541                   warninp_t wi)
3542 {
3543     char                      err_buf[256];
3544     int                       i, m, g, nmol, npct;
3545     gmx_bool                  bCharge, bAcc;
3546     real                      gdt_max, *mgrp, mt;
3547     rvec                      acc;
3548     gmx_mtop_atomloop_block_t aloopb;
3549     gmx_mtop_atomloop_all_t   aloop;
3550     t_atom                   *atom;
3551     ivec                      AbsRef;
3552     char                      warn_buf[STRLEN];
3553
3554     set_warning_line(wi, mdparin, -1);
3555
3556     if (EI_DYNAMICS(ir->eI) && !EI_SD(ir->eI) && ir->eI != eiBD &&
3557         ir->comm_mode == ecmNO &&
3558         !(absolute_reference(ir, sys, FALSE, AbsRef) || ir->nsteps <= 10))
3559     {
3560         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");
3561     }
3562
3563     /* Check for pressure coupling with absolute position restraints */
3564     if (ir->epc != epcNO && ir->refcoord_scaling == erscNO)
3565     {
3566         absolute_reference(ir, sys, TRUE, AbsRef);
3567         {
3568             for (m = 0; m < DIM; m++)
3569             {
3570                 if (AbsRef[m] && norm2(ir->compress[m]) > 0)
3571                 {
3572                     warning(wi, "You are using pressure coupling with absolute position restraints, this will give artifacts. Use the refcoord_scaling option.");
3573                     break;
3574                 }
3575             }
3576         }
3577     }
3578
3579     bCharge = FALSE;
3580     aloopb  = gmx_mtop_atomloop_block_init(sys);
3581     while (gmx_mtop_atomloop_block_next(aloopb, &atom, &nmol))
3582     {
3583         if (atom->q != 0 || atom->qB != 0)
3584         {
3585             bCharge = TRUE;
3586         }
3587     }
3588
3589     if (!bCharge)
3590     {
3591         if (EEL_FULL(ir->coulombtype))
3592         {
3593             sprintf(err_buf,
3594                     "You are using full electrostatics treatment %s for a system without charges.\n"
3595                     "This costs a lot of performance for just processing zeros, consider using %s instead.\n",
3596                     EELTYPE(ir->coulombtype), EELTYPE(eelCUT));
3597             warning(wi, err_buf);
3598         }
3599     }
3600     else
3601     {
3602         if (ir->coulombtype == eelCUT && ir->rcoulomb > 0 && !ir->implicit_solvent)
3603         {
3604             sprintf(err_buf,
3605                     "You are using a plain Coulomb cut-off, which might produce artifacts.\n"
3606                     "You might want to consider using %s electrostatics.\n",
3607                     EELTYPE(eelPME));
3608             warning_note(wi, err_buf);
3609         }
3610     }
3611
3612     /* Generalized reaction field */
3613     if (ir->opts.ngtc == 0)
3614     {
3615         sprintf(err_buf, "No temperature coupling while using coulombtype %s",
3616                 eel_names[eelGRF]);
3617         CHECK(ir->coulombtype == eelGRF);
3618     }
3619     else
3620     {
3621         sprintf(err_buf, "When using coulombtype = %s"
3622                 " ref-t for temperature coupling should be > 0",
3623                 eel_names[eelGRF]);
3624         CHECK((ir->coulombtype == eelGRF) && (ir->opts.ref_t[0] <= 0));
3625     }
3626
3627     if (ir->eI == eiSD1 &&
3628         (gmx_mtop_ftype_count(sys, F_CONSTR) > 0 ||
3629          gmx_mtop_ftype_count(sys, F_SETTLE) > 0))
3630     {
3631         sprintf(warn_buf, "With constraints integrator %s is less accurate, consider using %s instead", ei_names[ir->eI], ei_names[eiSD2]);
3632         warning_note(wi, warn_buf);
3633     }
3634
3635     bAcc = FALSE;
3636     for (i = 0; (i < sys->groups.grps[egcACC].nr); i++)
3637     {
3638         for (m = 0; (m < DIM); m++)
3639         {
3640             if (fabs(ir->opts.acc[i][m]) > 1e-6)
3641             {
3642                 bAcc = TRUE;
3643             }
3644         }
3645     }
3646     if (bAcc)
3647     {
3648         clear_rvec(acc);
3649         snew(mgrp, sys->groups.grps[egcACC].nr);
3650         aloop = gmx_mtop_atomloop_all_init(sys);
3651         while (gmx_mtop_atomloop_all_next(aloop, &i, &atom))
3652         {
3653             mgrp[ggrpnr(&sys->groups, egcACC, i)] += atom->m;
3654         }
3655         mt = 0.0;
3656         for (i = 0; (i < sys->groups.grps[egcACC].nr); i++)
3657         {
3658             for (m = 0; (m < DIM); m++)
3659             {
3660                 acc[m] += ir->opts.acc[i][m]*mgrp[i];
3661             }
3662             mt += mgrp[i];
3663         }
3664         for (m = 0; (m < DIM); m++)
3665         {
3666             if (fabs(acc[m]) > 1e-6)
3667             {
3668                 const char *dim[DIM] = { "X", "Y", "Z" };
3669                 fprintf(stderr,
3670                         "Net Acceleration in %s direction, will %s be corrected\n",
3671                         dim[m], ir->nstcomm != 0 ? "" : "not");
3672                 if (ir->nstcomm != 0 && m < ndof_com(ir))
3673                 {
3674                     acc[m] /= mt;
3675                     for (i = 0; (i < sys->groups.grps[egcACC].nr); i++)
3676                     {
3677                         ir->opts.acc[i][m] -= acc[m];
3678                     }
3679                 }
3680             }
3681         }
3682         sfree(mgrp);
3683     }
3684
3685     if (ir->efep != efepNO && ir->fepvals->sc_alpha != 0 &&
3686         !gmx_within_tol(sys->ffparams.reppow, 12.0, 10*GMX_DOUBLE_EPS))
3687     {
3688         gmx_fatal(FARGS, "Soft-core interactions are only supported with VdW repulsion power 12");
3689     }
3690
3691     if (ir->ePull != epullNO)
3692     {
3693         if (ir->pull->grp[0].nat == 0)
3694         {
3695             absolute_reference(ir, sys, FALSE, AbsRef);
3696             for (m = 0; m < DIM; m++)
3697             {
3698                 if (ir->pull->dim[m] && !AbsRef[m])
3699                 {
3700                     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.");
3701                     break;
3702                 }
3703             }
3704         }
3705
3706         if (ir->pull->eGeom == epullgDIRPBC)
3707         {
3708             for (i = 0; i < 3; i++)
3709             {
3710                 for (m = 0; m <= i; m++)
3711                 {
3712                     if ((ir->epc != epcNO && ir->compress[i][m] != 0) ||
3713                         ir->deform[i][m] != 0)
3714                     {
3715                         for (g = 1; g < ir->pull->ngrp; g++)
3716                         {
3717                             if (ir->pull->grp[g].vec[m] != 0)
3718                             {
3719                                 gmx_fatal(FARGS, "Can not have dynamic box while using pull geometry '%s' (dim %c)", EPULLGEOM(ir->pull->eGeom), 'x'+m);
3720                             }
3721                         }
3722                     }
3723                 }
3724             }
3725         }
3726     }
3727
3728     check_disre(sys);
3729 }
3730
3731 void double_check(t_inputrec *ir, matrix box, gmx_bool bConstr, warninp_t wi)
3732 {
3733     real        min_size;
3734     gmx_bool    bTWIN;
3735     char        warn_buf[STRLEN];
3736     const char *ptr;
3737
3738     ptr = check_box(ir->ePBC, box);
3739     if (ptr)
3740     {
3741         warning_error(wi, ptr);
3742     }
3743
3744     if (bConstr && ir->eConstrAlg == econtSHAKE)
3745     {
3746         if (ir->shake_tol <= 0.0)
3747         {
3748             sprintf(warn_buf, "ERROR: shake-tol must be > 0 instead of %g\n",
3749                     ir->shake_tol);
3750             warning_error(wi, warn_buf);
3751         }
3752
3753         if (IR_TWINRANGE(*ir) && ir->nstlist > 1)
3754         {
3755             sprintf(warn_buf, "With twin-range cut-off's and SHAKE the virial and the pressure are incorrect.");
3756             if (ir->epc == epcNO)
3757             {
3758                 warning(wi, warn_buf);
3759             }
3760             else
3761             {
3762                 warning_error(wi, warn_buf);
3763             }
3764         }
3765     }
3766
3767     if ( (ir->eConstrAlg == econtLINCS) && bConstr)
3768     {
3769         /* If we have Lincs constraints: */
3770         if (ir->eI == eiMD && ir->etc == etcNO &&
3771             ir->eConstrAlg == econtLINCS && ir->nLincsIter == 1)
3772         {
3773             sprintf(warn_buf, "For energy conservation with LINCS, lincs_iter should be 2 or larger.\n");
3774             warning_note(wi, warn_buf);
3775         }
3776
3777         if ((ir->eI == eiCG || ir->eI == eiLBFGS) && (ir->nProjOrder < 8))
3778         {
3779             sprintf(warn_buf, "For accurate %s with LINCS constraints, lincs-order should be 8 or more.", ei_names[ir->eI]);
3780             warning_note(wi, warn_buf);
3781         }
3782         if (ir->epc == epcMTTK)
3783         {
3784             warning_error(wi, "MTTK not compatible with lincs -- use shake instead.");
3785         }
3786     }
3787
3788     if (ir->LincsWarnAngle > 90.0)
3789     {
3790         sprintf(warn_buf, "lincs-warnangle can not be larger than 90 degrees, setting it to 90.\n");
3791         warning(wi, warn_buf);
3792         ir->LincsWarnAngle = 90.0;
3793     }
3794
3795     if (ir->ePBC != epbcNONE)
3796     {
3797         if (ir->nstlist == 0)
3798         {
3799             warning(wi, "With nstlist=0 atoms are only put into the box at step 0, therefore drifting atoms might cause the simulation to crash.");
3800         }
3801         bTWIN = (ir->rlistlong > ir->rlist);
3802         if (ir->ns_type == ensGRID)
3803         {
3804             if (sqr(ir->rlistlong) >= max_cutoff2(ir->ePBC, box))
3805             {
3806                 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",
3807                         bTWIN ? (ir->rcoulomb == ir->rlistlong ? "rcoulomb" : "rvdw") : "rlist");
3808                 warning_error(wi, warn_buf);
3809             }
3810         }
3811         else
3812         {
3813             min_size = min(box[XX][XX], min(box[YY][YY], box[ZZ][ZZ]));
3814             if (2*ir->rlistlong >= min_size)
3815             {
3816                 sprintf(warn_buf, "ERROR: One of the box lengths is smaller than twice the cut-off length. Increase the box size or decrease rlist.");
3817                 warning_error(wi, warn_buf);
3818                 if (TRICLINIC(box))
3819                 {
3820                     fprintf(stderr, "Grid search might allow larger cut-off's than simple search with triclinic boxes.");
3821                 }
3822             }
3823         }
3824     }
3825 }
3826
3827 void check_chargegroup_radii(const gmx_mtop_t *mtop, const t_inputrec *ir,
3828                              rvec *x,
3829                              warninp_t wi)
3830 {
3831     real rvdw1, rvdw2, rcoul1, rcoul2;
3832     char warn_buf[STRLEN];
3833
3834     calc_chargegroup_radii(mtop, x, &rvdw1, &rvdw2, &rcoul1, &rcoul2);
3835
3836     if (rvdw1 > 0)
3837     {
3838         printf("Largest charge group radii for Van der Waals: %5.3f, %5.3f nm\n",
3839                rvdw1, rvdw2);
3840     }
3841     if (rcoul1 > 0)
3842     {
3843         printf("Largest charge group radii for Coulomb:       %5.3f, %5.3f nm\n",
3844                rcoul1, rcoul2);
3845     }
3846
3847     if (ir->rlist > 0)
3848     {
3849         if (rvdw1  + rvdw2  > ir->rlist ||
3850             rcoul1 + rcoul2 > ir->rlist)
3851         {
3852             sprintf(warn_buf,
3853                     "The sum of the two largest charge group radii (%f) "
3854                     "is larger than rlist (%f)\n",
3855                     max(rvdw1+rvdw2, rcoul1+rcoul2), ir->rlist);
3856             warning(wi, warn_buf);
3857         }
3858         else
3859         {
3860             /* Here we do not use the zero at cut-off macro,
3861              * since user defined interactions might purposely
3862              * not be zero at the cut-off.
3863              */
3864             if ((EVDW_IS_ZERO_AT_CUTOFF(ir->vdwtype) ||
3865                  ir->vdw_modifier != eintmodNONE) &&
3866                 rvdw1 + rvdw2 > ir->rlistlong - ir->rvdw)
3867             {
3868                 sprintf(warn_buf, "The sum of the two largest charge group "
3869                         "radii (%f) is larger than %s (%f) - rvdw (%f).\n"
3870                         "With exact cut-offs, better performance can be "
3871                         "obtained with cutoff-scheme = %s, because it "
3872                         "does not use charge groups at all.",
3873                         rvdw1+rvdw2,
3874                         ir->rlistlong > ir->rlist ? "rlistlong" : "rlist",
3875                         ir->rlistlong, ir->rvdw,
3876                         ecutscheme_names[ecutsVERLET]);
3877                 if (ir_NVE(ir))
3878                 {
3879                     warning(wi, warn_buf);
3880                 }
3881                 else
3882                 {
3883                     warning_note(wi, warn_buf);
3884                 }
3885             }
3886             if ((EEL_IS_ZERO_AT_CUTOFF(ir->coulombtype) ||
3887                  ir->coulomb_modifier != eintmodNONE) &&
3888                 rcoul1 + rcoul2 > ir->rlistlong - ir->rcoulomb)
3889             {
3890                 sprintf(warn_buf, "The sum of the two largest charge group radii (%f) is larger than %s (%f) - rcoulomb (%f).\n"
3891                         "With exact cut-offs, better performance can be obtained with cutoff-scheme = %s, because it does not use charge groups at all.",
3892                         rcoul1+rcoul2,
3893                         ir->rlistlong > ir->rlist ? "rlistlong" : "rlist",
3894                         ir->rlistlong, ir->rcoulomb,
3895                         ecutscheme_names[ecutsVERLET]);
3896                 if (ir_NVE(ir))
3897                 {
3898                     warning(wi, warn_buf);
3899                 }
3900                 else
3901                 {
3902                     warning_note(wi, warn_buf);
3903                 }
3904             }
3905         }
3906     }
3907 }