src/gromacs/gmxpreprocess/readir.c

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