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