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             !((EEL_MIGHT_BE_ZERO_AT_CUTOFF(ir->coulombtype) && ir->rcoulomb > ir->rlist) ||
 278               (EVDW_MIGHT_BE_ZERO_AT_CUTOFF(ir->vdwtype)    && 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 (EEL_MIGHT_BE_ZERO_AT_CUTOFF(ir->coulombtype))
1080     {
1081         if (EEL_SWITCHED(ir->coulombtype))
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 (EVDW_MIGHT_BE_ZERO_AT_CUTOFF(ir->vdwtype))
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 (EVDW_SWITCHED(ir->vdwtype))
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 (EEL_IS_ZERO_AT_CUTOFF(ir->coulombtype)
1213             && (ir->rlistlong <= ir->rcoulomb))
1214         {
1215             sprintf(warn_buf, "For energy conservation with switch/shift potentials, %s should be 0.1 to 0.3 nm larger than rcoulomb.",
1216                     IR_TWINRANGE(*ir) ? "rlistlong" : "rlist");
1217             warning_note(wi, warn_buf);
1218         }
1219         if (EVDW_SWITCHED(ir->vdwtype) && (ir->rlistlong <= ir->rvdw))
1220         {
1221             sprintf(warn_buf, "For energy conservation with switch/shift potentials, %s should be 0.1 to 0.3 nm larger than rvdw.",
1222                     IR_TWINRANGE(*ir) ? "rlistlong" : "rlist");
1223             warning_note(wi, warn_buf);
1224         }
1225     }
1226
1227     if (ir->vdwtype == evdwUSER && ir->eDispCorr != edispcNO)
1228     {
1229         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.");
1230     }
1231
1232     if (ir->nstlist == -1)
1233     {
1234         sprintf(err_buf, "With nstlist=-1 rvdw and rcoulomb should be smaller than rlist to account for diffusion and possibly charge-group radii");
1235         CHECK(ir->rvdw >= ir->rlist || ir->rcoulomb >= ir->rlist);
1236     }
1237     sprintf(err_buf, "nstlist can not be smaller than -1");
1238     CHECK(ir->nstlist < -1);
1239
1240     if (ir->eI == eiLBFGS && (ir->coulombtype == eelCUT || ir->vdwtype == evdwCUT)
1241         && ir->rvdw != 0)
1242     {
1243         warning(wi, "For efficient BFGS minimization, use switch/shift/pme instead of cut-off.");
1244     }
1245
1246     if (ir->eI == eiLBFGS && ir->nbfgscorr <= 0)
1247     {
1248         warning(wi, "Using L-BFGS with nbfgscorr<=0 just gets you steepest descent.");
1249     }
1250
1251     /* ENERGY CONSERVATION */
1252     if (ir_NVE(ir) && ir->cutoff_scheme == ecutsGROUP)
1253     {
1254         if (!EVDW_MIGHT_BE_ZERO_AT_CUTOFF(ir->vdwtype) && ir->rvdw > 0 && ir->vdw_modifier == eintmodNONE)
1255         {
1256             sprintf(warn_buf, "You are using a cut-off for VdW interactions with NVE, for good energy conservation use vdwtype = %s (possibly with DispCorr)",
1257                     evdw_names[evdwSHIFT]);
1258             warning_note(wi, warn_buf);
1259         }
1260         if (!EEL_MIGHT_BE_ZERO_AT_CUTOFF(ir->coulombtype) && ir->rcoulomb > 0 && ir->coulomb_modifier == eintmodNONE)
1261         {
1262             sprintf(warn_buf, "You are using a cut-off for electrostatics with NVE, for good energy conservation use coulombtype = %s or %s",
1263                     eel_names[eelPMESWITCH], eel_names[eelRF_ZERO]);
1264             warning_note(wi, warn_buf);
1265         }
1266     }
1267
1268     /* IMPLICIT SOLVENT */
1269     if (ir->coulombtype == eelGB_NOTUSED)
1270     {
1271         ir->coulombtype      = eelCUT;
1272         ir->implicit_solvent = eisGBSA;
1273         fprintf(stderr, "Note: Old option for generalized born electrostatics given:\n"
1274                 "Changing coulombtype from \"generalized-born\" to \"cut-off\" and instead\n"
1275                 "setting implicit-solvent value to \"GBSA\" in input section.\n");
1276     }
1277
1278     if (ir->sa_algorithm == esaSTILL)
1279     {
1280         sprintf(err_buf, "Still SA algorithm not available yet, use %s or %s instead\n", esa_names[esaAPPROX], esa_names[esaNO]);
1281         CHECK(ir->sa_algorithm == esaSTILL);
1282     }
1283
1284     if (ir->implicit_solvent == eisGBSA)
1285     {
1286         sprintf(err_buf, "With GBSA implicit solvent, rgbradii must be equal to rlist.");
1287         CHECK(ir->rgbradii != ir->rlist);
1288
1289         if (ir->coulombtype != eelCUT)
1290         {
1291             sprintf(err_buf, "With GBSA, coulombtype must be equal to %s\n", eel_names[eelCUT]);
1292             CHECK(ir->coulombtype != eelCUT);
1293         }
1294         if (ir->vdwtype != evdwCUT)
1295         {
1296             sprintf(err_buf, "With GBSA, vdw-type must be equal to %s\n", evdw_names[evdwCUT]);
1297             CHECK(ir->vdwtype != evdwCUT);
1298         }
1299         if (ir->nstgbradii < 1)
1300         {
1301             sprintf(warn_buf, "Using GBSA with nstgbradii<1, setting nstgbradii=1");
1302             warning_note(wi, warn_buf);
1303             ir->nstgbradii = 1;
1304         }
1305         if (ir->sa_algorithm == esaNO)
1306         {
1307             sprintf(warn_buf, "No SA (non-polar) calculation requested together with GB. Are you sure this is what you want?\n");
1308             warning_note(wi, warn_buf);
1309         }
1310         if (ir->sa_surface_tension < 0 && ir->sa_algorithm != esaNO)
1311         {
1312             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");
1313             warning_note(wi, warn_buf);
1314
1315             if (ir->gb_algorithm == egbSTILL)
1316             {
1317                 ir->sa_surface_tension = 0.0049 * CAL2JOULE * 100;
1318             }
1319             else
1320             {
1321                 ir->sa_surface_tension = 0.0054 * CAL2JOULE * 100;
1322             }
1323         }
1324         if (ir->sa_surface_tension == 0 && ir->sa_algorithm != esaNO)
1325         {
1326             sprintf(err_buf, "Surface tension set to 0 while SA-calculation requested\n");
1327             CHECK(ir->sa_surface_tension == 0 && ir->sa_algorithm != esaNO);
1328         }
1329
1330     }
1331
1332     if (ir->bAdress)
1333     {
1334         if (ir->cutoff_scheme != ecutsGROUP)
1335         {
1336             warning_error(wi, "AdresS simulation supports only cutoff-scheme=group");
1337         }
1338         if (!EI_SD(ir->eI))
1339         {
1340             warning_error(wi, "AdresS simulation supports only stochastic dynamics");
1341         }
1342         if (ir->epc != epcNO)
1343         {
1344             warning_error(wi, "AdresS simulation does not support pressure coupling");
1345         }
1346         if (EEL_FULL(ir->coulombtype))
1347         {
1348             warning_error(wi, "AdresS simulation does not support long-range electrostatics");
1349         }
1350     }
1351 }
1352
1353 /* count the number of text elemets separated by whitespace in a string.
1354     str = the input string
1355     maxptr = the maximum number of allowed elements
1356     ptr = the output array of pointers to the first character of each element
1357     returns: the number of elements. */
1358 int str_nelem(const char *str, int maxptr, char *ptr[])
1359 {
1360     int   np = 0;
1361     char *copy0, *copy;
1362
1363     copy0 = strdup(str);
1364     copy  = copy0;
1365     ltrim(copy);
1366     while (*copy != '\0')
1367     {
1368         if (np >= maxptr)
1369         {
1370             gmx_fatal(FARGS, "Too many groups on line: '%s' (max is %d)",
1371                       str, maxptr);
1372         }
1373         if (ptr)
1374         {
1375             ptr[np] = copy;
1376         }
1377         np++;
1378         while ((*copy != '\0') && !isspace(*copy))
1379         {
1380             copy++;
1381         }
1382         if (*copy != '\0')
1383         {
1384             *copy = '\0';
1385             copy++;
1386         }
1387         ltrim(copy);
1388     }
1389     if (ptr == NULL)
1390     {
1391         sfree(copy0);
1392     }
1393
1394     return np;
1395 }
1396
1397 /* interpret a number of doubles from a string and put them in an array,
1398    after allocating space for them.
1399    str = the input string
1400    n = the (pre-allocated) number of doubles read
1401    r = the output array of doubles. */
1402 static void parse_n_real(char *str, int *n, real **r)
1403 {
1404     char *ptr[MAXPTR];
1405     int   i;
1406
1407     *n = str_nelem(str, MAXPTR, ptr);
1408
1409     snew(*r, *n);
1410     for (i = 0; i < *n; i++)
1411     {
1412         (*r)[i] = strtod(ptr[i], NULL);
1413     }
1414 }
1415
1416 static void do_fep_params(t_inputrec *ir, char fep_lambda[][STRLEN], char weights[STRLEN])
1417 {
1418
1419     int         i, j, max_n_lambda, nweights, nfep[efptNR];
1420     t_lambda   *fep    = ir->fepvals;
1421     t_expanded *expand = ir->expandedvals;
1422     real      **count_fep_lambdas;
1423     gmx_bool    bOneLambda = TRUE;
1424
1425     snew(count_fep_lambdas, efptNR);
1426
1427     /* FEP input processing */
1428     /* first, identify the number of lambda values for each type.
1429        All that are nonzero must have the same number */
1430
1431     for (i = 0; i < efptNR; i++)
1432     {
1433         parse_n_real(fep_lambda[i], &(nfep[i]), &(count_fep_lambdas[i]));
1434     }
1435
1436     /* now, determine the number of components.  All must be either zero, or equal. */
1437
1438     max_n_lambda = 0;
1439     for (i = 0; i < efptNR; i++)
1440     {
1441         if (nfep[i] > max_n_lambda)
1442         {
1443             max_n_lambda = nfep[i];  /* here's a nonzero one.  All of them
1444                                         must have the same number if its not zero.*/
1445             break;
1446         }
1447     }
1448
1449     for (i = 0; i < efptNR; i++)
1450     {
1451         if (nfep[i] == 0)
1452         {
1453             ir->fepvals->separate_dvdl[i] = FALSE;
1454         }
1455         else if (nfep[i] == max_n_lambda)
1456         {
1457             if (i != efptTEMPERATURE)  /* we treat this differently -- not really a reason to compute the derivative with
1458                                           respect to the temperature currently */
1459             {
1460                 ir->fepvals->separate_dvdl[i] = TRUE;
1461             }
1462         }
1463         else
1464         {
1465             gmx_fatal(FARGS, "Number of lambdas (%d) for FEP type %s not equal to number of other types (%d)",
1466                       nfep[i], efpt_names[i], max_n_lambda);
1467         }
1468     }
1469     /* we don't print out dhdl if the temperature is changing, since we can't correctly define dhdl in this case */
1470     ir->fepvals->separate_dvdl[efptTEMPERATURE] = FALSE;
1471
1472     /* the number of lambdas is the number we've read in, which is either zero
1473        or the same for all */
1474     fep->n_lambda = max_n_lambda;
1475
1476     /* allocate space for the array of lambda values */
1477     snew(fep->all_lambda, efptNR);
1478     /* if init_lambda is defined, we need to set lambda */
1479     if ((fep->init_lambda > 0) && (fep->n_lambda == 0))
1480     {
1481         ir->fepvals->separate_dvdl[efptFEP] = TRUE;
1482     }
1483     /* otherwise allocate the space for all of the lambdas, and transfer the data */
1484     for (i = 0; i < efptNR; i++)
1485     {
1486         snew(fep->all_lambda[i], fep->n_lambda);
1487         if (nfep[i] > 0)  /* if it's zero, then the count_fep_lambda arrays
1488                              are zero */
1489         {
1490             for (j = 0; j < fep->n_lambda; j++)
1491             {
1492                 fep->all_lambda[i][j] = (double)count_fep_lambdas[i][j];
1493             }
1494             sfree(count_fep_lambdas[i]);
1495         }
1496     }
1497     sfree(count_fep_lambdas);
1498
1499     /* "fep-vals" is either zero or the full number. If zero, we'll need to define fep-lambdas for internal
1500        bookkeeping -- for now, init_lambda */
1501
1502     if ((nfep[efptFEP] == 0) && (fep->init_lambda >= 0))
1503     {
1504         for (i = 0; i < fep->n_lambda; i++)
1505         {
1506             fep->all_lambda[efptFEP][i] = fep->init_lambda;
1507         }
1508     }
1509
1510     /* check to see if only a single component lambda is defined, and soft core is defined.
1511        In this case, turn on coulomb soft core */
1512
1513     if (max_n_lambda == 0)
1514     {
1515         bOneLambda = TRUE;
1516     }
1517     else
1518     {
1519         for (i = 0; i < efptNR; i++)
1520         {
1521             if ((nfep[i] != 0) && (i != efptFEP))
1522             {
1523                 bOneLambda = FALSE;
1524             }
1525         }
1526     }
1527     if ((bOneLambda) && (fep->sc_alpha > 0))
1528     {
1529         fep->bScCoul = TRUE;
1530     }
1531
1532     /* Fill in the others with the efptFEP if they are not explicitly
1533        specified (i.e. nfep[i] == 0).  This means if fep is not defined,
1534        they are all zero. */
1535
1536     for (i = 0; i < efptNR; i++)
1537     {
1538         if ((nfep[i] == 0) && (i != efptFEP))
1539         {
1540             for (j = 0; j < fep->n_lambda; j++)
1541             {
1542                 fep->all_lambda[i][j] = fep->all_lambda[efptFEP][j];
1543             }
1544         }
1545     }
1546
1547
1548     /* make it easier if sc_r_power = 48 by increasing it to the 4th power, to be in the right scale. */
1549     if (fep->sc_r_power == 48)
1550     {
1551         if (fep->sc_alpha > 0.1)
1552         {
1553             gmx_fatal(FARGS, "sc_alpha (%f) for sc_r_power = 48 should usually be between 0.001 and 0.004", fep->sc_alpha);
1554         }
1555     }
1556
1557     expand = ir->expandedvals;
1558     /* now read in the weights */
1559     parse_n_real(weights, &nweights, &(expand->init_lambda_weights));
1560     if (nweights == 0)
1561     {
1562         snew(expand->init_lambda_weights, fep->n_lambda); /* initialize to zero */
1563     }
1564     else if (nweights != fep->n_lambda)
1565     {
1566         gmx_fatal(FARGS, "Number of weights (%d) is not equal to number of lambda values (%d)",
1567                   nweights, fep->n_lambda);
1568     }
1569     if ((expand->nstexpanded < 0) && (ir->efep != efepNO))
1570     {
1571         expand->nstexpanded = fep->nstdhdl;
1572         /* if you don't specify nstexpanded when doing expanded ensemble free energy calcs, it is set to nstdhdl */
1573     }
1574     if ((expand->nstexpanded < 0) && ir->bSimTemp)
1575     {
1576         expand->nstexpanded = 2*(int)(ir->opts.tau_t[0]/ir->delta_t);
1577         /* if you don't specify nstexpanded when doing expanded ensemble simulated tempering, it is set to
1578            2*tau_t just to be careful so it's not to frequent  */
1579     }
1580 }
1581
1582
1583 static void do_simtemp_params(t_inputrec *ir)
1584 {
1585
1586     snew(ir->simtempvals->temperatures, ir->fepvals->n_lambda);
1587     GetSimTemps(ir->fepvals->n_lambda, ir->simtempvals, ir->fepvals->all_lambda[efptTEMPERATURE]);
1588
1589     return;
1590 }
1591
1592 static void do_wall_params(t_inputrec *ir,
1593                            char *wall_atomtype, char *wall_density,
1594                            t_gromppopts *opts)
1595 {
1596     int    nstr, i;
1597     char  *names[MAXPTR];
1598     double dbl;
1599
1600     opts->wall_atomtype[0] = NULL;
1601     opts->wall_atomtype[1] = NULL;
1602
1603     ir->wall_atomtype[0] = -1;
1604     ir->wall_atomtype[1] = -1;
1605     ir->wall_density[0]  = 0;
1606     ir->wall_density[1]  = 0;
1607
1608     if (ir->nwall > 0)
1609     {
1610         nstr = str_nelem(wall_atomtype, MAXPTR, names);
1611         if (nstr != ir->nwall)
1612         {
1613             gmx_fatal(FARGS, "Expected %d elements for wall_atomtype, found %d",
1614                       ir->nwall, nstr);
1615         }
1616         for (i = 0; i < ir->nwall; i++)
1617         {
1618             opts->wall_atomtype[i] = strdup(names[i]);
1619         }
1620
1621         if (ir->wall_type == ewt93 || ir->wall_type == ewt104)
1622         {
1623             nstr = str_nelem(wall_density, MAXPTR, names);
1624             if (nstr != ir->nwall)
1625             {
1626                 gmx_fatal(FARGS, "Expected %d elements for wall-density, found %d", ir->nwall, nstr);
1627             }
1628             for (i = 0; i < ir->nwall; i++)
1629             {
1630                 sscanf(names[i], "%lf", &dbl);
1631                 if (dbl <= 0)
1632                 {
1633                     gmx_fatal(FARGS, "wall-density[%d] = %f\n", i, dbl);
1634                 }
1635                 ir->wall_density[i] = dbl;
1636             }
1637         }
1638     }
1639 }
1640
1641 static void add_wall_energrps(gmx_groups_t *groups, int nwall, t_symtab *symtab)
1642 {
1643     int     i;
1644     t_grps *grps;
1645     char    str[STRLEN];
1646
1647     if (nwall > 0)
1648     {
1649         srenew(groups->grpname, groups->ngrpname+nwall);
1650         grps = &(groups->grps[egcENER]);
1651         srenew(grps->nm_ind, grps->nr+nwall);
1652         for (i = 0; i < nwall; i++)
1653         {
1654             sprintf(str, "wall%d", i);
1655             groups->grpname[groups->ngrpname] = put_symtab(symtab, str);
1656             grps->nm_ind[grps->nr++]          = groups->ngrpname++;
1657         }
1658     }
1659 }
1660
1661 void read_expandedparams(int *ninp_p, t_inpfile **inp_p,
1662                          t_expanded *expand, warninp_t wi)
1663 {
1664     int        ninp, nerror = 0;
1665     t_inpfile *inp;
1666
1667     ninp   = *ninp_p;
1668     inp    = *inp_p;
1669
1670     /* read expanded ensemble parameters */
1671     CCTYPE ("expanded ensemble variables");
1672     ITYPE ("nstexpanded", expand->nstexpanded, -1);
1673     EETYPE("lmc-stats", expand->elamstats, elamstats_names);
1674     EETYPE("lmc-move", expand->elmcmove, elmcmove_names);
1675     EETYPE("lmc-weights-equil", expand->elmceq, elmceq_names);
1676     ITYPE ("weight-equil-number-all-lambda", expand->equil_n_at_lam, -1);
1677     ITYPE ("weight-equil-number-samples", expand->equil_samples, -1);
1678     ITYPE ("weight-equil-number-steps", expand->equil_steps, -1);
1679     RTYPE ("weight-equil-wl-delta", expand->equil_wl_delta, -1);
1680     RTYPE ("weight-equil-count-ratio", expand->equil_ratio, -1);
1681     CCTYPE("Seed for Monte Carlo in lambda space");
1682     ITYPE ("lmc-seed", expand->lmc_seed, -1);
1683     RTYPE ("mc-temperature", expand->mc_temp, -1);
1684     ITYPE ("lmc-repeats", expand->lmc_repeats, 1);
1685     ITYPE ("lmc-gibbsdelta", expand->gibbsdeltalam, -1);
1686     ITYPE ("lmc-forced-nstart", expand->lmc_forced_nstart, 0);
1687     EETYPE("symmetrized-transition-matrix", expand->bSymmetrizedTMatrix, yesno_names);
1688     ITYPE("nst-transition-matrix", expand->nstTij, -1);
1689     ITYPE ("mininum-var-min", expand->minvarmin, 100); /*default is reasonable */
1690     ITYPE ("weight-c-range", expand->c_range, 0);      /* default is just C=0 */
1691     RTYPE ("wl-scale", expand->wl_scale, 0.8);
1692     RTYPE ("wl-ratio", expand->wl_ratio, 0.8);
1693     RTYPE ("init-wl-delta", expand->init_wl_delta, 1.0);
1694     EETYPE("wl-oneovert", expand->bWLoneovert, yesno_names);
1695
1696     *ninp_p   = ninp;
1697     *inp_p    = inp;
1698
1699     return;
1700 }
1701
1702 void get_ir(const char *mdparin, const char *mdparout,
1703             t_inputrec *ir, t_gromppopts *opts,
1704             warninp_t wi)
1705 {
1706     char       *dumstr[2];
1707     double      dumdub[2][6];
1708     t_inpfile  *inp;
1709     const char *tmp;
1710     int         i, j, m, ninp;
1711     char        warn_buf[STRLEN];
1712     t_lambda   *fep    = ir->fepvals;
1713     t_expanded *expand = ir->expandedvals;
1714
1715     init_inputrec_strings();
1716     inp = read_inpfile(mdparin, &ninp, wi);
1717
1718     snew(dumstr[0], STRLEN);
1719     snew(dumstr[1], STRLEN);
1720
1721     if (-1 == search_einp(ninp, inp, "cutoff-scheme"))
1722     {
1723         sprintf(warn_buf,
1724                 "%s did not specify a value for the .mdp option "
1725                 "\"cutoff-scheme\". Probably it was first intended for use "
1726                 "with GROMACS before 4.6. In 4.6, the Verlet scheme was "
1727                 "introduced, but the group scheme was still the default. "
1728                 "The default is now the Verlet scheme, so you will observe "
1729                 "different behaviour.", mdparin);
1730         warning_note(wi, warn_buf);
1731     }
1732
1733     /* remove the following deprecated commands */
1734     REM_TYPE("title");
1735     REM_TYPE("cpp");
1736     REM_TYPE("domain-decomposition");
1737     REM_TYPE("andersen-seed");
1738     REM_TYPE("dihre");
1739     REM_TYPE("dihre-fc");
1740     REM_TYPE("dihre-tau");
1741     REM_TYPE("nstdihreout");
1742     REM_TYPE("nstcheckpoint");
1743
1744     /* replace the following commands with the clearer new versions*/
1745     REPL_TYPE("unconstrained-start", "continuation");
1746     REPL_TYPE("foreign-lambda", "fep-lambdas");
1747     REPL_TYPE("verlet-buffer-drift", "verlet-buffer-tolerance");
1748     REPL_TYPE("nstxtcout", "nstxout-compressed");
1749     REPL_TYPE("xtc-grps", "compressed-x-grps");
1750     REPL_TYPE("xtc-precision", "compressed-x-precision");
1751
1752     CCTYPE ("VARIOUS PREPROCESSING OPTIONS");
1753     CTYPE ("Preprocessor information: use cpp syntax.");
1754     CTYPE ("e.g.: -I/home/joe/doe -I/home/mary/roe");
1755     STYPE ("include", opts->include,  NULL);
1756     CTYPE ("e.g.: -DPOSRES -DFLEXIBLE (note these variable names are case sensitive)");
1757     STYPE ("define",  opts->define,   NULL);
1758
1759     CCTYPE ("RUN CONTROL PARAMETERS");
1760     EETYPE("integrator",  ir->eI,         ei_names);
1761     CTYPE ("Start time and timestep in ps");
1762     RTYPE ("tinit",   ir->init_t, 0.0);
1763     RTYPE ("dt",      ir->delta_t,    0.001);
1764     STEPTYPE ("nsteps",   ir->nsteps,     0);
1765     CTYPE ("For exact run continuation or redoing part of a run");
1766     STEPTYPE ("init-step", ir->init_step,  0);
1767     CTYPE ("Part index is updated automatically on checkpointing (keeps files separate)");
1768     ITYPE ("simulation-part", ir->simulation_part, 1);
1769     CTYPE ("mode for center of mass motion removal");
1770     EETYPE("comm-mode",   ir->comm_mode,  ecm_names);
1771     CTYPE ("number of steps for center of mass motion removal");
1772     ITYPE ("nstcomm", ir->nstcomm,    100);
1773     CTYPE ("group(s) for center of mass motion removal");
1774     STYPE ("comm-grps",   is->vcm,            NULL);
1775
1776     CCTYPE ("LANGEVIN DYNAMICS OPTIONS");
1777     CTYPE ("Friction coefficient (amu/ps) and random seed");
1778     RTYPE ("bd-fric",     ir->bd_fric,    0.0);
1779     ITYPE ("ld-seed",     ir->ld_seed,    -1);
1780
1781     /* Em stuff */
1782     CCTYPE ("ENERGY MINIMIZATION OPTIONS");
1783     CTYPE ("Force tolerance and initial step-size");
1784     RTYPE ("emtol",       ir->em_tol,     10.0);
1785     RTYPE ("emstep",      ir->em_stepsize, 0.01);
1786     CTYPE ("Max number of iterations in relax-shells");
1787     ITYPE ("niter",       ir->niter,      20);
1788     CTYPE ("Step size (ps^2) for minimization of flexible constraints");
1789     RTYPE ("fcstep",      ir->fc_stepsize, 0);
1790     CTYPE ("Frequency of steepest descents steps when doing CG");
1791     ITYPE ("nstcgsteep",  ir->nstcgsteep, 1000);
1792     ITYPE ("nbfgscorr",   ir->nbfgscorr,  10);
1793
1794     CCTYPE ("TEST PARTICLE INSERTION OPTIONS");
1795     RTYPE ("rtpi",    ir->rtpi,   0.05);
1796
1797     /* Output options */
1798     CCTYPE ("OUTPUT CONTROL OPTIONS");
1799     CTYPE ("Output frequency for coords (x), velocities (v) and forces (f)");
1800     ITYPE ("nstxout", ir->nstxout,    0);
1801     ITYPE ("nstvout", ir->nstvout,    0);
1802     ITYPE ("nstfout", ir->nstfout,    0);
1803     ir->nstcheckpoint = 1000;
1804     CTYPE ("Output frequency for energies to log file and energy file");
1805     ITYPE ("nstlog",  ir->nstlog, 1000);
1806     ITYPE ("nstcalcenergy", ir->nstcalcenergy, 100);
1807     ITYPE ("nstenergy",   ir->nstenergy,  1000);
1808     CTYPE ("Output frequency and precision for .xtc file");
1809     ITYPE ("nstxout-compressed", ir->nstxout_compressed,  0);
1810     RTYPE ("compressed-x-precision", ir->x_compression_precision, 1000.0);
1811     CTYPE ("This selects the subset of atoms for the compressed");
1812     CTYPE ("trajectory file. You can select multiple groups. By");
1813     CTYPE ("default, all atoms will be written.");
1814     STYPE ("compressed-x-grps", is->x_compressed_groups, NULL);
1815     CTYPE ("Selection of energy groups");
1816     STYPE ("energygrps",  is->energy,         NULL);
1817
1818     /* Neighbor searching */
1819     CCTYPE ("NEIGHBORSEARCHING PARAMETERS");
1820     CTYPE ("cut-off scheme (Verlet: particle based cut-offs, group: using charge groups)");
1821     EETYPE("cutoff-scheme",     ir->cutoff_scheme,    ecutscheme_names);
1822     CTYPE ("nblist update frequency");
1823     ITYPE ("nstlist", ir->nstlist,    10);
1824     CTYPE ("ns algorithm (simple or grid)");
1825     EETYPE("ns-type",     ir->ns_type,    ens_names);
1826     /* set ndelta to the optimal value of 2 */
1827     ir->ndelta = 2;
1828     CTYPE ("Periodic boundary conditions: xyz, no, xy");
1829     EETYPE("pbc",         ir->ePBC,       epbc_names);
1830     EETYPE("periodic-molecules", ir->bPeriodicMols, yesno_names);
1831     CTYPE ("Allowed energy error due to the Verlet buffer in kJ/mol/ps per atom,");
1832     CTYPE ("a value of -1 means: use rlist");
1833     RTYPE("verlet-buffer-tolerance", ir->verletbuf_tol,    0.005);
1834     CTYPE ("nblist cut-off");
1835     RTYPE ("rlist",   ir->rlist,  1.0);
1836     CTYPE ("long-range cut-off for switched potentials");
1837     RTYPE ("rlistlong",   ir->rlistlong,  -1);
1838     ITYPE ("nstcalclr",   ir->nstcalclr,  -1);
1839
1840     /* Electrostatics */
1841     CCTYPE ("OPTIONS FOR ELECTROSTATICS AND VDW");
1842     CTYPE ("Method for doing electrostatics");
1843     EETYPE("coulombtype", ir->coulombtype,    eel_names);
1844     EETYPE("coulomb-modifier",    ir->coulomb_modifier,    eintmod_names);
1845     CTYPE ("cut-off lengths");
1846     RTYPE ("rcoulomb-switch", ir->rcoulomb_switch,    0.0);
1847     RTYPE ("rcoulomb",    ir->rcoulomb,   1.0);
1848     CTYPE ("Relative dielectric constant for the medium and the reaction field");
1849     RTYPE ("epsilon-r",   ir->epsilon_r,  1.0);
1850     RTYPE ("epsilon-rf",  ir->epsilon_rf, 0.0);
1851     CTYPE ("Method for doing Van der Waals");
1852     EETYPE("vdw-type",    ir->vdwtype,    evdw_names);
1853     EETYPE("vdw-modifier",    ir->vdw_modifier,    eintmod_names);
1854     CTYPE ("cut-off lengths");
1855     RTYPE ("rvdw-switch", ir->rvdw_switch,    0.0);
1856     RTYPE ("rvdw",    ir->rvdw,   1.0);
1857     CTYPE ("Apply long range dispersion corrections for Energy and Pressure");
1858     EETYPE("DispCorr",    ir->eDispCorr,  edispc_names);
1859     CTYPE ("Extension of the potential lookup tables beyond the cut-off");
1860     RTYPE ("table-extension", ir->tabext, 1.0);
1861     CTYPE ("Separate tables between energy group pairs");
1862     STYPE ("energygrp-table", is->egptable,   NULL);
1863     CTYPE ("Spacing for the PME/PPPM FFT grid");
1864     RTYPE ("fourierspacing", ir->fourier_spacing, 0.12);
1865     CTYPE ("FFT grid size, when a value is 0 fourierspacing will be used");
1866     ITYPE ("fourier-nx",  ir->nkx,         0);
1867     ITYPE ("fourier-ny",  ir->nky,         0);
1868     ITYPE ("fourier-nz",  ir->nkz,         0);
1869     CTYPE ("EWALD/PME/PPPM parameters");
1870     ITYPE ("pme-order",   ir->pme_order,   4);
1871     RTYPE ("ewald-rtol",  ir->ewald_rtol, 0.00001);
1872     RTYPE ("ewald-rtol-lj", ir->ewald_rtol_lj, 0.001);
1873     EETYPE("lj-pme-comb-rule", ir->ljpme_combination_rule, eljpme_names);
1874     EETYPE("ewald-geometry", ir->ewald_geometry, eewg_names);
1875     RTYPE ("epsilon-surface", ir->epsilon_surface, 0.0);
1876     EETYPE("optimize-fft", ir->bOptFFT,  yesno_names);
1877
1878     CCTYPE("IMPLICIT SOLVENT ALGORITHM");
1879     EETYPE("implicit-solvent", ir->implicit_solvent, eis_names);
1880
1881     CCTYPE ("GENERALIZED BORN ELECTROSTATICS");
1882     CTYPE ("Algorithm for calculating Born radii");
1883     EETYPE("gb-algorithm", ir->gb_algorithm, egb_names);
1884     CTYPE ("Frequency of calculating the Born radii inside rlist");
1885     ITYPE ("nstgbradii", ir->nstgbradii, 1);
1886     CTYPE ("Cutoff for Born radii calculation; the contribution from atoms");
1887     CTYPE ("between rlist and rgbradii is updated every nstlist steps");
1888     RTYPE ("rgbradii",  ir->rgbradii, 1.0);
1889     CTYPE ("Dielectric coefficient of the implicit solvent");
1890     RTYPE ("gb-epsilon-solvent", ir->gb_epsilon_solvent, 80.0);
1891     CTYPE ("Salt concentration in M for Generalized Born models");
1892     RTYPE ("gb-saltconc",  ir->gb_saltconc, 0.0);
1893     CTYPE ("Scaling factors used in the OBC GB model. Default values are OBC(II)");
1894     RTYPE ("gb-obc-alpha", ir->gb_obc_alpha, 1.0);
1895     RTYPE ("gb-obc-beta", ir->gb_obc_beta, 0.8);
1896     RTYPE ("gb-obc-gamma", ir->gb_obc_gamma, 4.85);
1897     RTYPE ("gb-dielectric-offset", ir->gb_dielectric_offset, 0.009);
1898     EETYPE("sa-algorithm", ir->sa_algorithm, esa_names);
1899     CTYPE ("Surface tension (kJ/mol/nm^2) for the SA (nonpolar surface) part of GBSA");
1900     CTYPE ("The value -1 will set default value for Still/HCT/OBC GB-models.");
1901     RTYPE ("sa-surface-tension", ir->sa_surface_tension, -1);
1902
1903     /* Coupling stuff */
1904     CCTYPE ("OPTIONS FOR WEAK COUPLING ALGORITHMS");
1905     CTYPE ("Temperature coupling");
1906     EETYPE("tcoupl",  ir->etc,        etcoupl_names);
1907     ITYPE ("nsttcouple", ir->nsttcouple,  -1);
1908     ITYPE("nh-chain-length",     ir->opts.nhchainlength, 10);
1909     EETYPE("print-nose-hoover-chain-variables", ir->bPrintNHChains, yesno_names);
1910     CTYPE ("Groups to couple separately");
1911     STYPE ("tc-grps",     is->tcgrps,         NULL);
1912     CTYPE ("Time constant (ps) and reference temperature (K)");
1913     STYPE ("tau-t",   is->tau_t,      NULL);
1914     STYPE ("ref-t",   is->ref_t,      NULL);
1915     CTYPE ("pressure coupling");
1916     EETYPE("pcoupl",  ir->epc,        epcoupl_names);
1917     EETYPE("pcoupltype",  ir->epct,       epcoupltype_names);
1918     ITYPE ("nstpcouple", ir->nstpcouple,  -1);
1919     CTYPE ("Time constant (ps), compressibility (1/bar) and reference P (bar)");
1920     RTYPE ("tau-p",   ir->tau_p,  1.0);
1921     STYPE ("compressibility", dumstr[0],  NULL);
1922     STYPE ("ref-p",       dumstr[1],      NULL);
1923     CTYPE ("Scaling of reference coordinates, No, All or COM");
1924     EETYPE ("refcoord-scaling", ir->refcoord_scaling, erefscaling_names);
1925
1926     /* QMMM */
1927     CCTYPE ("OPTIONS FOR QMMM calculations");
1928     EETYPE("QMMM", ir->bQMMM, yesno_names);
1929     CTYPE ("Groups treated Quantum Mechanically");
1930     STYPE ("QMMM-grps",  is->QMMM,          NULL);
1931     CTYPE ("QM method");
1932     STYPE("QMmethod",     is->QMmethod, NULL);
1933     CTYPE ("QMMM scheme");
1934     EETYPE("QMMMscheme",  ir->QMMMscheme,    eQMMMscheme_names);
1935     CTYPE ("QM basisset");
1936     STYPE("QMbasis",      is->QMbasis, NULL);
1937     CTYPE ("QM charge");
1938     STYPE ("QMcharge",    is->QMcharge, NULL);
1939     CTYPE ("QM multiplicity");
1940     STYPE ("QMmult",      is->QMmult, NULL);
1941     CTYPE ("Surface Hopping");
1942     STYPE ("SH",          is->bSH, NULL);
1943     CTYPE ("CAS space options");
1944     STYPE ("CASorbitals",      is->CASorbitals,   NULL);
1945     STYPE ("CASelectrons",     is->CASelectrons,  NULL);
1946     STYPE ("SAon", is->SAon, NULL);
1947     STYPE ("SAoff", is->SAoff, NULL);
1948     STYPE ("SAsteps", is->SAsteps, NULL);
1949     CTYPE ("Scale factor for MM charges");
1950     RTYPE ("MMChargeScaleFactor", ir->scalefactor, 1.0);
1951     CTYPE ("Optimization of QM subsystem");
1952     STYPE ("bOPT",          is->bOPT, NULL);
1953     STYPE ("bTS",          is->bTS, NULL);
1954
1955     /* Simulated annealing */
1956     CCTYPE("SIMULATED ANNEALING");
1957     CTYPE ("Type of annealing for each temperature group (no/single/periodic)");
1958     STYPE ("annealing",   is->anneal,      NULL);
1959     CTYPE ("Number of time points to use for specifying annealing in each group");
1960     STYPE ("annealing-npoints", is->anneal_npoints, NULL);
1961     CTYPE ("List of times at the annealing points for each group");
1962     STYPE ("annealing-time",       is->anneal_time,       NULL);
1963     CTYPE ("Temp. at each annealing point, for each group.");
1964     STYPE ("annealing-temp",  is->anneal_temp,  NULL);
1965
1966     /* Startup run */
1967     CCTYPE ("GENERATE VELOCITIES FOR STARTUP RUN");
1968     EETYPE("gen-vel",     opts->bGenVel,  yesno_names);
1969     RTYPE ("gen-temp",    opts->tempi,    300.0);
1970     ITYPE ("gen-seed",    opts->seed,     -1);
1971
1972     /* Shake stuff */
1973     CCTYPE ("OPTIONS FOR BONDS");
1974     EETYPE("constraints", opts->nshake,   constraints);
1975     CTYPE ("Type of constraint algorithm");
1976     EETYPE("constraint-algorithm",  ir->eConstrAlg, econstr_names);
1977     CTYPE ("Do not constrain the start configuration");
1978     EETYPE("continuation", ir->bContinuation, yesno_names);
1979     CTYPE ("Use successive overrelaxation to reduce the number of shake iterations");
1980     EETYPE("Shake-SOR", ir->bShakeSOR, yesno_names);
1981     CTYPE ("Relative tolerance of shake");
1982     RTYPE ("shake-tol", ir->shake_tol, 0.0001);
1983     CTYPE ("Highest order in the expansion of the constraint coupling matrix");
1984     ITYPE ("lincs-order", ir->nProjOrder, 4);
1985     CTYPE ("Number of iterations in the final step of LINCS. 1 is fine for");
1986     CTYPE ("normal simulations, but use 2 to conserve energy in NVE runs.");
1987     CTYPE ("For energy minimization with constraints it should be 4 to 8.");
1988     ITYPE ("lincs-iter", ir->nLincsIter, 1);
1989     CTYPE ("Lincs will write a warning to the stderr if in one step a bond");
1990     CTYPE ("rotates over more degrees than");
1991     RTYPE ("lincs-warnangle", ir->LincsWarnAngle, 30.0);
1992     CTYPE ("Convert harmonic bonds to morse potentials");
1993     EETYPE("morse",       opts->bMorse, yesno_names);
1994
1995     /* Energy group exclusions */
1996     CCTYPE ("ENERGY GROUP EXCLUSIONS");
1997     CTYPE ("Pairs of energy groups for which all non-bonded interactions are excluded");
1998     STYPE ("energygrp-excl", is->egpexcl,     NULL);
1999
2000     /* Walls */
2001     CCTYPE ("WALLS");
2002     CTYPE ("Number of walls, type, atom types, densities and box-z scale factor for Ewald");
2003     ITYPE ("nwall", ir->nwall, 0);
2004     EETYPE("wall-type",     ir->wall_type,   ewt_names);
2005     RTYPE ("wall-r-linpot", ir->wall_r_linpot, -1);
2006     STYPE ("wall-atomtype", is->wall_atomtype, NULL);
2007     STYPE ("wall-density",  is->wall_density,  NULL);
2008     RTYPE ("wall-ewald-zfac", ir->wall_ewald_zfac, 3);
2009
2010     /* COM pulling */
2011     CCTYPE("COM PULLING");
2012     CTYPE("Pull type: no, umbrella, constraint or constant-force");
2013     EETYPE("pull",          ir->ePull, epull_names);
2014     if (ir->ePull != epullNO)
2015     {
2016         snew(ir->pull, 1);
2017         is->pull_grp = read_pullparams(&ninp, &inp, ir->pull, &opts->pull_start, wi);
2018     }
2019
2020     /* Enforced rotation */
2021     CCTYPE("ENFORCED ROTATION");
2022     CTYPE("Enforced rotation: No or Yes");
2023     EETYPE("rotation",       ir->bRot, yesno_names);
2024     if (ir->bRot)
2025     {
2026         snew(ir->rot, 1);
2027         is->rot_grp = read_rotparams(&ninp, &inp, ir->rot, wi);
2028     }
2029
2030     /* Refinement */
2031     CCTYPE("NMR refinement stuff");
2032     CTYPE ("Distance restraints type: No, Simple or Ensemble");
2033     EETYPE("disre",       ir->eDisre,     edisre_names);
2034     CTYPE ("Force weighting of pairs in one distance restraint: Conservative or Equal");
2035     EETYPE("disre-weighting", ir->eDisreWeighting, edisreweighting_names);
2036     CTYPE ("Use sqrt of the time averaged times the instantaneous violation");
2037     EETYPE("disre-mixed", ir->bDisreMixed, yesno_names);
2038     RTYPE ("disre-fc",    ir->dr_fc,  1000.0);
2039     RTYPE ("disre-tau",   ir->dr_tau, 0.0);
2040     CTYPE ("Output frequency for pair distances to energy file");
2041     ITYPE ("nstdisreout", ir->nstdisreout, 100);
2042     CTYPE ("Orientation restraints: No or Yes");
2043     EETYPE("orire",       opts->bOrire,   yesno_names);
2044     CTYPE ("Orientation restraints force constant and tau for time averaging");
2045     RTYPE ("orire-fc",    ir->orires_fc,  0.0);
2046     RTYPE ("orire-tau",   ir->orires_tau, 0.0);
2047     STYPE ("orire-fitgrp", is->orirefitgrp,    NULL);
2048     CTYPE ("Output frequency for trace(SD) and S to energy file");
2049     ITYPE ("nstorireout", ir->nstorireout, 100);
2050
2051     /* free energy variables */
2052     CCTYPE ("Free energy variables");
2053     EETYPE("free-energy", ir->efep, efep_names);
2054     STYPE ("couple-moltype",  is->couple_moltype,  NULL);
2055     EETYPE("couple-lambda0", opts->couple_lam0, couple_lam);
2056     EETYPE("couple-lambda1", opts->couple_lam1, couple_lam);
2057     EETYPE("couple-intramol", opts->bCoupleIntra, yesno_names);
2058
2059     RTYPE ("init-lambda", fep->init_lambda, -1); /* start with -1 so
2060                                                     we can recognize if
2061                                                     it was not entered */
2062     ITYPE ("init-lambda-state", fep->init_fep_state, -1);
2063     RTYPE ("delta-lambda", fep->delta_lambda, 0.0);
2064     ITYPE ("nstdhdl", fep->nstdhdl, 50);
2065     STYPE ("fep-lambdas", is->fep_lambda[efptFEP], NULL);
2066     STYPE ("mass-lambdas", is->fep_lambda[efptMASS], NULL);
2067     STYPE ("coul-lambdas", is->fep_lambda[efptCOUL], NULL);
2068     STYPE ("vdw-lambdas", is->fep_lambda[efptVDW], NULL);
2069     STYPE ("bonded-lambdas", is->fep_lambda[efptBONDED], NULL);
2070     STYPE ("restraint-lambdas", is->fep_lambda[efptRESTRAINT], NULL);
2071     STYPE ("temperature-lambdas", is->fep_lambda[efptTEMPERATURE], NULL);
2072     ITYPE ("calc-lambda-neighbors", fep->lambda_neighbors, 1);
2073     STYPE ("init-lambda-weights", is->lambda_weights, NULL);
2074     EETYPE("dhdl-print-energy", fep->bPrintEnergy, yesno_names);
2075     RTYPE ("sc-alpha", fep->sc_alpha, 0.0);
2076     ITYPE ("sc-power", fep->sc_power, 1);
2077     RTYPE ("sc-r-power", fep->sc_r_power, 6.0);
2078     RTYPE ("sc-sigma", fep->sc_sigma, 0.3);
2079     EETYPE("sc-coul", fep->bScCoul, yesno_names);
2080     ITYPE ("dh_hist_size", fep->dh_hist_size, 0);
2081     RTYPE ("dh_hist_spacing", fep->dh_hist_spacing, 0.1);
2082     EETYPE("separate-dhdl-file", fep->separate_dhdl_file,
2083            separate_dhdl_file_names);
2084     EETYPE("dhdl-derivatives", fep->dhdl_derivatives, dhdl_derivatives_names);
2085     ITYPE ("dh_hist_size", fep->dh_hist_size, 0);
2086     RTYPE ("dh_hist_spacing", fep->dh_hist_spacing, 0.1);
2087
2088     /* Non-equilibrium MD stuff */
2089     CCTYPE("Non-equilibrium MD stuff");
2090     STYPE ("acc-grps",    is->accgrps,        NULL);
2091     STYPE ("accelerate",  is->acc,            NULL);
2092     STYPE ("freezegrps",  is->freeze,         NULL);
2093     STYPE ("freezedim",   is->frdim,          NULL);
2094     RTYPE ("cos-acceleration", ir->cos_accel, 0);
2095     STYPE ("deform",      is->deform,         NULL);
2096
2097     /* simulated tempering variables */
2098     CCTYPE("simulated tempering variables");
2099     EETYPE("simulated-tempering", ir->bSimTemp, yesno_names);
2100     EETYPE("simulated-tempering-scaling", ir->simtempvals->eSimTempScale, esimtemp_names);
2101     RTYPE("sim-temp-low", ir->simtempvals->simtemp_low, 300.0);
2102     RTYPE("sim-temp-high", ir->simtempvals->simtemp_high, 300.0);
2103
2104     /* expanded ensemble variables */
2105     if (ir->efep == efepEXPANDED || ir->bSimTemp)
2106     {
2107         read_expandedparams(&ninp, &inp, expand, wi);
2108     }
2109
2110     /* Electric fields */
2111     CCTYPE("Electric fields");
2112     CTYPE ("Format is number of terms (int) and for all terms an amplitude (real)");
2113     CTYPE ("and a phase angle (real)");
2114     STYPE ("E-x",     is->efield_x,   NULL);
2115     STYPE ("E-xt",    is->efield_xt,  NULL);
2116     STYPE ("E-y",     is->efield_y,   NULL);
2117     STYPE ("E-yt",    is->efield_yt,  NULL);
2118     STYPE ("E-z",     is->efield_z,   NULL);
2119     STYPE ("E-zt",    is->efield_zt,  NULL);
2120
2121     CCTYPE("Ion/water position swapping for computational electrophysiology setups");
2122     CTYPE("Swap positions along direction: no, X, Y, Z");
2123     EETYPE("swapcoords", ir->eSwapCoords, eSwapTypes_names);
2124     if (ir->eSwapCoords != eswapNO)
2125     {
2126         snew(ir->swap, 1);
2127         CTYPE("Swap attempt frequency");
2128         ITYPE("swap-frequency", ir->swap->nstswap, 1);
2129         CTYPE("Two index groups that contain the compartment-partitioning atoms");
2130         STYPE("split-group0", splitgrp0, NULL);
2131         STYPE("split-group1", splitgrp1, NULL);
2132         CTYPE("Use center of mass of split groups (yes/no), otherwise center of geometry is used");
2133         EETYPE("massw-split0", ir->swap->massw_split[0], yesno_names);
2134         EETYPE("massw-split1", ir->swap->massw_split[1], yesno_names);
2135
2136         CTYPE("Group name of ions that can be exchanged with solvent molecules");
2137         STYPE("swap-group", swapgrp, NULL);
2138         CTYPE("Group name of solvent molecules");
2139         STYPE("solvent-group", solgrp, NULL);
2140
2141         CTYPE("Split cylinder: radius, upper and lower extension (nm) (this will define the channels)");
2142         CTYPE("Note that the split cylinder settings do not have an influence on the swapping protocol,");
2143         CTYPE("however, if correctly defined, the ion permeation events are counted per channel");
2144         RTYPE("cyl0-r", ir->swap->cyl0r, 2.0);
2145         RTYPE("cyl0-up", ir->swap->cyl0u, 1.0);
2146         RTYPE("cyl0-down", ir->swap->cyl0l, 1.0);
2147         RTYPE("cyl1-r", ir->swap->cyl1r, 2.0);
2148         RTYPE("cyl1-up", ir->swap->cyl1u, 1.0);
2149         RTYPE("cyl1-down", ir->swap->cyl1l, 1.0);
2150
2151         CTYPE("Average the number of ions per compartment over these many swap attempt steps");
2152         ITYPE("coupl-steps", ir->swap->nAverage, 10);
2153         CTYPE("Requested number of anions and cations for each of the two compartments");
2154         CTYPE("-1 means fix the numbers as found in time step 0");
2155         ITYPE("anionsA", ir->swap->nanions[0], -1);
2156         ITYPE("cationsA", ir->swap->ncations[0], -1);
2157         ITYPE("anionsB", ir->swap->nanions[1], -1);
2158         ITYPE("cationsB", ir->swap->ncations[1], -1);
2159         CTYPE("Start to swap ions if threshold difference to requested count is reached");
2160         RTYPE("threshold", ir->swap->threshold, 1.0);
2161     }
2162
2163     /* AdResS defined thingies */
2164     CCTYPE ("AdResS parameters");
2165     EETYPE("adress",       ir->bAdress, yesno_names);
2166     if (ir->bAdress)
2167     {
2168         snew(ir->adress, 1);
2169         read_adressparams(&ninp, &inp, ir->adress, wi);
2170     }
2171
2172     /* User defined thingies */
2173     CCTYPE ("User defined thingies");
2174     STYPE ("user1-grps",  is->user1,          NULL);
2175     STYPE ("user2-grps",  is->user2,          NULL);
2176     ITYPE ("userint1",    ir->userint1,   0);
2177     ITYPE ("userint2",    ir->userint2,   0);
2178     ITYPE ("userint3",    ir->userint3,   0);
2179     ITYPE ("userint4",    ir->userint4,   0);
2180     RTYPE ("userreal1",   ir->userreal1,  0);
2181     RTYPE ("userreal2",   ir->userreal2,  0);
2182     RTYPE ("userreal3",   ir->userreal3,  0);
2183     RTYPE ("userreal4",   ir->userreal4,  0);
2184 #undef CTYPE
2185
2186     write_inpfile(mdparout, ninp, inp, FALSE, wi);
2187     for (i = 0; (i < ninp); i++)
2188     {
2189         sfree(inp[i].name);
2190         sfree(inp[i].value);
2191     }
2192     sfree(inp);
2193
2194     /* Process options if necessary */
2195     for (m = 0; m < 2; m++)
2196     {
2197         for (i = 0; i < 2*DIM; i++)
2198         {
2199             dumdub[m][i] = 0.0;
2200         }
2201         if (ir->epc)
2202         {
2203             switch (ir->epct)
2204             {
2205                 case epctISOTROPIC:
2206                     if (sscanf(dumstr[m], "%lf", &(dumdub[m][XX])) != 1)
2207                     {
2208                         warning_error(wi, "Pressure coupling not enough values (I need 1)");
2209                     }
2210                     dumdub[m][YY] = dumdub[m][ZZ] = dumdub[m][XX];
2211                     break;
2212                 case epctSEMIISOTROPIC:
2213                 case epctSURFACETENSION:
2214                     if (sscanf(dumstr[m], "%lf%lf",
2215                                &(dumdub[m][XX]), &(dumdub[m][ZZ])) != 2)
2216                     {
2217                         warning_error(wi, "Pressure coupling not enough values (I need 2)");
2218                     }
2219                     dumdub[m][YY] = dumdub[m][XX];
2220                     break;
2221                 case epctANISOTROPIC:
2222                     if (sscanf(dumstr[m], "%lf%lf%lf%lf%lf%lf",
2223                                &(dumdub[m][XX]), &(dumdub[m][YY]), &(dumdub[m][ZZ]),
2224                                &(dumdub[m][3]), &(dumdub[m][4]), &(dumdub[m][5])) != 6)
2225                     {
2226                         warning_error(wi, "Pressure coupling not enough values (I need 6)");
2227                     }
2228                     break;
2229                 default:
2230                     gmx_fatal(FARGS, "Pressure coupling type %s not implemented yet",
2231                               epcoupltype_names[ir->epct]);
2232             }
2233         }
2234     }
2235     clear_mat(ir->ref_p);
2236     clear_mat(ir->compress);
2237     for (i = 0; i < DIM; i++)
2238     {
2239         ir->ref_p[i][i]    = dumdub[1][i];
2240         ir->compress[i][i] = dumdub[0][i];
2241     }
2242     if (ir->epct == epctANISOTROPIC)
2243     {
2244         ir->ref_p[XX][YY] = dumdub[1][3];
2245         ir->ref_p[XX][ZZ] = dumdub[1][4];
2246         ir->ref_p[YY][ZZ] = dumdub[1][5];
2247         if (ir->ref_p[XX][YY] != 0 && ir->ref_p[XX][ZZ] != 0 && ir->ref_p[YY][ZZ] != 0)
2248         {
2249             warning(wi, "All off-diagonal reference pressures are non-zero. Are you sure you want to apply a threefold shear stress?\n");
2250         }
2251         ir->compress[XX][YY] = dumdub[0][3];
2252         ir->compress[XX][ZZ] = dumdub[0][4];
2253         ir->compress[YY][ZZ] = dumdub[0][5];
2254         for (i = 0; i < DIM; i++)
2255         {
2256             for (m = 0; m < i; m++)
2257             {
2258                 ir->ref_p[i][m]    = ir->ref_p[m][i];
2259                 ir->compress[i][m] = ir->compress[m][i];
2260             }
2261         }
2262     }
2263
2264     if (ir->comm_mode == ecmNO)
2265     {
2266         ir->nstcomm = 0;
2267     }
2268
2269     opts->couple_moltype = NULL;
2270     if (strlen(is->couple_moltype) > 0)
2271     {
2272         if (ir->efep != efepNO)
2273         {
2274             opts->couple_moltype = strdup(is->couple_moltype);
2275             if (opts->couple_lam0 == opts->couple_lam1)
2276             {
2277                 warning(wi, "The lambda=0 and lambda=1 states for coupling are identical");
2278             }
2279             if (ir->eI == eiMD && (opts->couple_lam0 == ecouplamNONE ||
2280                                    opts->couple_lam1 == ecouplamNONE))
2281             {
2282                 warning(wi, "For proper sampling of the (nearly) decoupled state, stochastic dynamics should be used");
2283             }
2284         }
2285         else
2286         {
2287             warning(wi, "Can not couple a molecule with free_energy = no");
2288         }
2289     }
2290     /* FREE ENERGY AND EXPANDED ENSEMBLE OPTIONS */
2291     if (ir->efep != efepNO)
2292     {
2293         if (fep->delta_lambda > 0)
2294         {
2295             ir->efep = efepSLOWGROWTH;
2296         }
2297     }
2298
2299     if (ir->bSimTemp)
2300     {
2301         fep->bPrintEnergy = TRUE;
2302         /* always print out the energy to dhdl if we are doing expanded ensemble, since we need the total energy
2303            if the temperature is changing. */
2304     }
2305
2306     if ((ir->efep != efepNO) || ir->bSimTemp)
2307     {
2308         ir->bExpanded = FALSE;
2309         if ((ir->efep == efepEXPANDED) || ir->bSimTemp)
2310         {
2311             ir->bExpanded = TRUE;
2312         }
2313         do_fep_params(ir, is->fep_lambda, is->lambda_weights);
2314         if (ir->bSimTemp) /* done after fep params */
2315         {
2316             do_simtemp_params(ir);
2317         }
2318     }
2319     else
2320     {
2321         ir->fepvals->n_lambda = 0;
2322     }
2323
2324     /* WALL PARAMETERS */
2325
2326     do_wall_params(ir, is->wall_atomtype, is->wall_density, opts);
2327
2328     /* ORIENTATION RESTRAINT PARAMETERS */
2329
2330     if (opts->bOrire && str_nelem(is->orirefitgrp, MAXPTR, NULL) != 1)
2331     {
2332         warning_error(wi, "ERROR: Need one orientation restraint fit group\n");
2333     }
2334
2335     /* DEFORMATION PARAMETERS */
2336
2337     clear_mat(ir->deform);
2338     for (i = 0; i < 6; i++)
2339     {
2340         dumdub[0][i] = 0;
2341     }
2342     m = sscanf(is->deform, "%lf %lf %lf %lf %lf %lf",
2343                &(dumdub[0][0]), &(dumdub[0][1]), &(dumdub[0][2]),
2344                &(dumdub[0][3]), &(dumdub[0][4]), &(dumdub[0][5]));
2345     for (i = 0; i < 3; i++)
2346     {
2347         ir->deform[i][i] = dumdub[0][i];
2348     }
2349     ir->deform[YY][XX] = dumdub[0][3];
2350     ir->deform[ZZ][XX] = dumdub[0][4];
2351     ir->deform[ZZ][YY] = dumdub[0][5];
2352     if (ir->epc != epcNO)
2353     {
2354         for (i = 0; i < 3; i++)
2355         {
2356             for (j = 0; j <= i; j++)
2357             {
2358                 if (ir->deform[i][j] != 0 && ir->compress[i][j] != 0)
2359                 {
2360                     warning_error(wi, "A box element has deform set and compressibility > 0");
2361                 }
2362             }
2363         }
2364         for (i = 0; i < 3; i++)
2365         {
2366             for (j = 0; j < i; j++)
2367             {
2368                 if (ir->deform[i][j] != 0)
2369                 {
2370                     for (m = j; m < DIM; m++)
2371                     {
2372                         if (ir->compress[m][j] != 0)
2373                         {
2374                             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.");
2375                             warning(wi, warn_buf);
2376                         }
2377                     }
2378                 }
2379             }
2380         }
2381     }
2382
2383     /* Ion/water position swapping checks */
2384     if (ir->eSwapCoords != eswapNO)
2385     {
2386         if (ir->swap->nstswap < 1)
2387         {
2388             warning_error(wi, "swap_frequency must be 1 or larger when ion swapping is requested");
2389         }
2390         if (ir->swap->nAverage < 1)
2391         {
2392             warning_error(wi, "coupl_steps must be 1 or larger.\n");
2393         }
2394         if (ir->swap->threshold < 1.0)
2395         {
2396             warning_error(wi, "Ion count threshold must be at least 1.\n");
2397         }
2398     }
2399
2400     sfree(dumstr[0]);
2401     sfree(dumstr[1]);
2402 }
2403
2404 static int search_QMstring(const char *s, int ng, const char *gn[])
2405 {
2406     /* same as normal search_string, but this one searches QM strings */
2407     int i;
2408
2409     for (i = 0; (i < ng); i++)
2410     {
2411         if (gmx_strcasecmp(s, gn[i]) == 0)
2412         {
2413             return i;
2414         }
2415     }
2416
2417     gmx_fatal(FARGS, "this QM method or basisset (%s) is not implemented\n!", s);
2418
2419     return -1;
2420
2421 } /* search_QMstring */
2422
2423 /* We would like gn to be const as well, but C doesn't allow this */
2424 int search_string(const char *s, int ng, char *gn[])
2425 {
2426     int i;
2427
2428     for (i = 0; (i < ng); i++)
2429     {
2430         if (gmx_strcasecmp(s, gn[i]) == 0)
2431         {
2432             return i;
2433         }
2434     }
2435
2436     gmx_fatal(FARGS,
2437               "Group %s referenced in the .mdp file was not found in the index file.\n"
2438               "Group names must match either [moleculetype] names or custom index group\n"
2439               "names, in which case you must supply an index file to the '-n' option\n"
2440               "of grompp.",
2441               s);
2442
2443     return -1;
2444 }
2445
2446 static gmx_bool do_numbering(int natoms, gmx_groups_t *groups, int ng, char *ptrs[],
2447                              t_blocka *block, char *gnames[],
2448                              int gtype, int restnm,
2449                              int grptp, gmx_bool bVerbose,
2450                              warninp_t wi)
2451 {
2452     unsigned short *cbuf;
2453     t_grps         *grps = &(groups->grps[gtype]);
2454     int             i, j, gid, aj, ognr, ntot = 0;
2455     const char     *title;
2456     gmx_bool        bRest;
2457     char            warn_buf[STRLEN];
2458
2459     if (debug)
2460     {
2461         fprintf(debug, "Starting numbering %d groups of type %d\n", ng, gtype);
2462     }
2463
2464     title = gtypes[gtype];
2465
2466     snew(cbuf, natoms);
2467     /* Mark all id's as not set */
2468     for (i = 0; (i < natoms); i++)
2469     {
2470         cbuf[i] = NOGID;
2471     }
2472
2473     snew(grps->nm_ind, ng+1); /* +1 for possible rest group */
2474     for (i = 0; (i < ng); i++)
2475     {
2476         /* Lookup the group name in the block structure */
2477         gid = search_string(ptrs[i], block->nr, gnames);
2478         if ((grptp != egrptpONE) || (i == 0))
2479         {
2480             grps->nm_ind[grps->nr++] = gid;
2481         }
2482         if (debug)
2483         {
2484             fprintf(debug, "Found gid %d for group %s\n", gid, ptrs[i]);
2485         }
2486
2487         /* Now go over the atoms in the group */
2488         for (j = block->index[gid]; (j < block->index[gid+1]); j++)
2489         {
2490
2491             aj = block->a[j];
2492
2493             /* Range checking */
2494             if ((aj < 0) || (aj >= natoms))
2495             {
2496                 gmx_fatal(FARGS, "Invalid atom number %d in indexfile", aj);
2497             }
2498             /* Lookup up the old group number */
2499             ognr = cbuf[aj];
2500             if (ognr != NOGID)
2501             {
2502                 gmx_fatal(FARGS, "Atom %d in multiple %s groups (%d and %d)",
2503                           aj+1, title, ognr+1, i+1);
2504             }
2505             else
2506             {
2507                 /* Store the group number in buffer */
2508                 if (grptp == egrptpONE)
2509                 {
2510                     cbuf[aj] = 0;
2511                 }
2512                 else
2513                 {
2514                     cbuf[aj] = i;
2515                 }
2516                 ntot++;
2517             }
2518         }
2519     }
2520
2521     /* Now check whether we have done all atoms */
2522     bRest = FALSE;
2523     if (ntot != natoms)
2524     {
2525         if (grptp == egrptpALL)
2526         {
2527             gmx_fatal(FARGS, "%d atoms are not part of any of the %s groups",
2528                       natoms-ntot, title);
2529         }
2530         else if (grptp == egrptpPART)
2531         {
2532             sprintf(warn_buf, "%d atoms are not part of any of the %s groups",
2533                     natoms-ntot, title);
2534             warning_note(wi, warn_buf);
2535         }
2536         /* Assign all atoms currently unassigned to a rest group */
2537         for (j = 0; (j < natoms); j++)
2538         {
2539             if (cbuf[j] == NOGID)
2540             {
2541                 cbuf[j] = grps->nr;
2542                 bRest   = TRUE;
2543             }
2544         }
2545         if (grptp != egrptpPART)
2546         {
2547             if (bVerbose)
2548             {
2549                 fprintf(stderr,
2550                         "Making dummy/rest group for %s containing %d elements\n",
2551                         title, natoms-ntot);
2552             }
2553             /* Add group name "rest" */
2554             grps->nm_ind[grps->nr] = restnm;
2555
2556             /* Assign the rest name to all atoms not currently assigned to a group */
2557             for (j = 0; (j < natoms); j++)
2558             {
2559                 if (cbuf[j] == NOGID)
2560                 {
2561                     cbuf[j] = grps->nr;
2562                 }
2563             }
2564             grps->nr++;
2565         }
2566     }
2567
2568     if (grps->nr == 1 && (ntot == 0 || ntot == natoms))
2569     {
2570         /* All atoms are part of one (or no) group, no index required */
2571         groups->ngrpnr[gtype] = 0;
2572         groups->grpnr[gtype]  = NULL;
2573     }
2574     else
2575     {
2576         groups->ngrpnr[gtype] = natoms;
2577         snew(groups->grpnr[gtype], natoms);
2578         for (j = 0; (j < natoms); j++)
2579         {
2580             groups->grpnr[gtype][j] = cbuf[j];
2581         }
2582     }
2583
2584     sfree(cbuf);
2585
2586     return (bRest && grptp == egrptpPART);
2587 }
2588
2589 static void calc_nrdf(gmx_mtop_t *mtop, t_inputrec *ir, char **gnames)
2590 {
2591     t_grpopts              *opts;
2592     gmx_groups_t           *groups;
2593     t_pull                 *pull;
2594     int                     natoms, ai, aj, i, j, d, g, imin, jmin;
2595     t_iatom                *ia;
2596     int                    *nrdf2, *na_vcm, na_tot;
2597     double                 *nrdf_tc, *nrdf_vcm, nrdf_uc, n_sub = 0;
2598     gmx_mtop_atomloop_all_t aloop;
2599     t_atom                 *atom;
2600     int                     mb, mol, ftype, as;
2601     gmx_molblock_t         *molb;
2602     gmx_moltype_t          *molt;
2603
2604     /* Calculate nrdf.
2605      * First calc 3xnr-atoms for each group
2606      * then subtract half a degree of freedom for each constraint
2607      *
2608      * Only atoms and nuclei contribute to the degrees of freedom...
2609      */
2610
2611     opts = &ir->opts;
2612
2613     groups = &mtop->groups;
2614     natoms = mtop->natoms;
2615
2616     /* Allocate one more for a possible rest group */
2617     /* We need to sum degrees of freedom into doubles,
2618      * since floats give too low nrdf's above 3 million atoms.
2619      */
2620     snew(nrdf_tc, groups->grps[egcTC].nr+1);
2621     snew(nrdf_vcm, groups->grps[egcVCM].nr+1);
2622     snew(na_vcm, groups->grps[egcVCM].nr+1);
2623
2624     for (i = 0; i < groups->grps[egcTC].nr; i++)
2625     {
2626         nrdf_tc[i] = 0;
2627     }
2628     for (i = 0; i < groups->grps[egcVCM].nr+1; i++)
2629     {
2630         nrdf_vcm[i] = 0;
2631     }
2632
2633     snew(nrdf2, natoms);
2634     aloop = gmx_mtop_atomloop_all_init(mtop);
2635     while (gmx_mtop_atomloop_all_next(aloop, &i, &atom))
2636     {
2637         nrdf2[i] = 0;
2638         if (atom->ptype == eptAtom || atom->ptype == eptNucleus)
2639         {
2640             g = ggrpnr(groups, egcFREEZE, i);
2641             /* Double count nrdf for particle i */
2642             for (d = 0; d < DIM; d++)
2643             {
2644                 if (opts->nFreeze[g][d] == 0)
2645                 {
2646                     nrdf2[i] += 2;
2647                 }
2648             }
2649             nrdf_tc [ggrpnr(groups, egcTC, i)]  += 0.5*nrdf2[i];
2650             nrdf_vcm[ggrpnr(groups, egcVCM, i)] += 0.5*nrdf2[i];
2651         }
2652     }
2653
2654     as = 0;
2655     for (mb = 0; mb < mtop->nmolblock; mb++)
2656     {
2657         molb = &mtop->molblock[mb];
2658         molt = &mtop->moltype[molb->type];
2659         atom = molt->atoms.atom;
2660         for (mol = 0; mol < molb->nmol; mol++)
2661         {
2662             for (ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++)
2663             {
2664                 ia = molt->ilist[ftype].iatoms;
2665                 for (i = 0; i < molt->ilist[ftype].nr; )
2666                 {
2667                     /* Subtract degrees of freedom for the constraints,
2668                      * if the particles still have degrees of freedom left.
2669                      * If one of the particles is a vsite or a shell, then all
2670                      * constraint motion will go there, but since they do not
2671                      * contribute to the constraints the degrees of freedom do not
2672                      * change.
2673                      */
2674                     ai = as + ia[1];
2675                     aj = as + ia[2];
2676                     if (((atom[ia[1]].ptype == eptNucleus) ||
2677                          (atom[ia[1]].ptype == eptAtom)) &&
2678                         ((atom[ia[2]].ptype == eptNucleus) ||
2679                          (atom[ia[2]].ptype == eptAtom)))
2680                     {
2681                         if (nrdf2[ai] > 0)
2682                         {
2683                             jmin = 1;
2684                         }
2685                         else
2686                         {
2687                             jmin = 2;
2688                         }
2689                         if (nrdf2[aj] > 0)
2690                         {
2691                             imin = 1;
2692                         }
2693                         else
2694                         {
2695                             imin = 2;
2696                         }
2697                         imin       = min(imin, nrdf2[ai]);
2698                         jmin       = min(jmin, nrdf2[aj]);
2699                         nrdf2[ai] -= imin;
2700                         nrdf2[aj] -= jmin;
2701                         nrdf_tc [ggrpnr(groups, egcTC, ai)]  -= 0.5*imin;
2702                         nrdf_tc [ggrpnr(groups, egcTC, aj)]  -= 0.5*jmin;
2703                         nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5*imin;
2704                         nrdf_vcm[ggrpnr(groups, egcVCM, aj)] -= 0.5*jmin;
2705                     }
2706                     ia += interaction_function[ftype].nratoms+1;
2707                     i  += interaction_function[ftype].nratoms+1;
2708                 }
2709             }
2710             ia = molt->ilist[F_SETTLE].iatoms;
2711             for (i = 0; i < molt->ilist[F_SETTLE].nr; )
2712             {
2713                 /* Subtract 1 dof from every atom in the SETTLE */
2714                 for (j = 0; j < 3; j++)
2715                 {
2716                     ai         = as + ia[1+j];
2717                     imin       = min(2, nrdf2[ai]);
2718                     nrdf2[ai] -= imin;
2719                     nrdf_tc [ggrpnr(groups, egcTC, ai)]  -= 0.5*imin;
2720                     nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5*imin;
2721                 }
2722                 ia += 4;
2723                 i  += 4;
2724             }
2725             as += molt->atoms.nr;
2726         }
2727     }
2728
2729     if (ir->ePull == epullCONSTRAINT)
2730     {
2731         /* Correct nrdf for the COM constraints.
2732          * We correct using the TC and VCM group of the first atom
2733          * in the reference and pull group. If atoms in one pull group
2734          * belong to different TC or VCM groups it is anyhow difficult
2735          * to determine the optimal nrdf assignment.
2736          */
2737         pull = ir->pull;
2738
2739         for (i = 0; i < pull->ncoord; i++)
2740         {
2741             imin = 1;
2742
2743             for (j = 0; j < 2; j++)
2744             {
2745                 const t_pull_group *pgrp;
2746
2747                 pgrp = &pull->group[pull->coord[i].group[j]];
2748
2749                 if (pgrp->nat > 0)
2750                 {
2751                     /* Subtract 1/2 dof from each group */
2752                     ai = pgrp->ind[0];
2753                     nrdf_tc [ggrpnr(groups, egcTC, ai)]  -= 0.5*imin;
2754                     nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5*imin;
2755                     if (nrdf_tc[ggrpnr(groups, egcTC, ai)] < 0)
2756                     {
2757                         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)]]);
2758                     }
2759                 }
2760                 else
2761                 {
2762                     /* We need to subtract the whole DOF from group j=1 */
2763                     imin += 1;
2764                 }
2765             }
2766         }
2767     }
2768
2769     if (ir->nstcomm != 0)
2770     {
2771         /* Subtract 3 from the number of degrees of freedom in each vcm group
2772          * when com translation is removed and 6 when rotation is removed
2773          * as well.
2774          */
2775         switch (ir->comm_mode)
2776         {
2777             case ecmLINEAR:
2778                 n_sub = ndof_com(ir);
2779                 break;
2780             case ecmANGULAR:
2781                 n_sub = 6;
2782                 break;
2783             default:
2784                 n_sub = 0;
2785                 gmx_incons("Checking comm_mode");
2786         }
2787
2788         for (i = 0; i < groups->grps[egcTC].nr; i++)
2789         {
2790             /* Count the number of atoms of TC group i for every VCM group */
2791             for (j = 0; j < groups->grps[egcVCM].nr+1; j++)
2792             {
2793                 na_vcm[j] = 0;
2794             }
2795             na_tot = 0;
2796             for (ai = 0; ai < natoms; ai++)
2797             {
2798                 if (ggrpnr(groups, egcTC, ai) == i)
2799                 {
2800                     na_vcm[ggrpnr(groups, egcVCM, ai)]++;
2801                     na_tot++;
2802                 }
2803             }
2804             /* Correct for VCM removal according to the fraction of each VCM
2805              * group present in this TC group.
2806              */
2807             nrdf_uc = nrdf_tc[i];
2808             if (debug)
2809             {
2810                 fprintf(debug, "T-group[%d] nrdf_uc = %g, n_sub = %g\n",
2811                         i, nrdf_uc, n_sub);
2812             }
2813             nrdf_tc[i] = 0;
2814             for (j = 0; j < groups->grps[egcVCM].nr+1; j++)
2815             {
2816                 if (nrdf_vcm[j] > n_sub)
2817                 {
2818                     nrdf_tc[i] += nrdf_uc*((double)na_vcm[j]/(double)na_tot)*
2819                         (nrdf_vcm[j] - n_sub)/nrdf_vcm[j];
2820                 }
2821                 if (debug)
2822                 {
2823                     fprintf(debug, "  nrdf_vcm[%d] = %g, nrdf = %g\n",
2824                             j, nrdf_vcm[j], nrdf_tc[i]);
2825                 }
2826             }
2827         }
2828     }
2829     for (i = 0; (i < groups->grps[egcTC].nr); i++)
2830     {
2831         opts->nrdf[i] = nrdf_tc[i];
2832         if (opts->nrdf[i] < 0)
2833         {
2834             opts->nrdf[i] = 0;
2835         }
2836         fprintf(stderr,
2837                 "Number of degrees of freedom in T-Coupling group %s is %.2f\n",
2838                 gnames[groups->grps[egcTC].nm_ind[i]], opts->nrdf[i]);
2839     }
2840
2841     sfree(nrdf2);
2842     sfree(nrdf_tc);
2843     sfree(nrdf_vcm);
2844     sfree(na_vcm);
2845 }
2846
2847 static void decode_cos(char *s, t_cosines *cosine)
2848 {
2849     char   *t;
2850     char    format[STRLEN], f1[STRLEN];
2851     double  a, phi;
2852     int     i;
2853
2854     t = strdup(s);
2855     trim(t);
2856
2857     cosine->n   = 0;
2858     cosine->a   = NULL;
2859     cosine->phi = NULL;
2860     if (strlen(t))
2861     {
2862         sscanf(t, "%d", &(cosine->n));
2863         if (cosine->n <= 0)
2864         {
2865             cosine->n = 0;
2866         }
2867         else
2868         {
2869             snew(cosine->a, cosine->n);
2870             snew(cosine->phi, cosine->n);
2871
2872             sprintf(format, "%%*d");
2873             for (i = 0; (i < cosine->n); i++)
2874             {
2875                 strcpy(f1, format);
2876                 strcat(f1, "%lf%lf");
2877                 if (sscanf(t, f1, &a, &phi) < 2)
2878                 {
2879                     gmx_fatal(FARGS, "Invalid input for electric field shift: '%s'", t);
2880                 }
2881                 cosine->a[i]   = a;
2882                 cosine->phi[i] = phi;
2883                 strcat(format, "%*lf%*lf");
2884             }
2885         }
2886     }
2887     sfree(t);
2888 }
2889
2890 static gmx_bool do_egp_flag(t_inputrec *ir, gmx_groups_t *groups,
2891                             const char *option, const char *val, int flag)
2892 {
2893     /* The maximum number of energy group pairs would be MAXPTR*(MAXPTR+1)/2.
2894      * But since this is much larger than STRLEN, such a line can not be parsed.
2895      * The real maximum is the number of names that fit in a string: STRLEN/2.
2896      */
2897 #define EGP_MAX (STRLEN/2)
2898     int      nelem, i, j, k, nr;
2899     char    *names[EGP_MAX];
2900     char  ***gnames;
2901     gmx_bool bSet;
2902
2903     gnames = groups->grpname;
2904
2905     nelem = str_nelem(val, EGP_MAX, names);
2906     if (nelem % 2 != 0)
2907     {
2908         gmx_fatal(FARGS, "The number of groups for %s is odd", option);
2909     }
2910     nr   = groups->grps[egcENER].nr;
2911     bSet = FALSE;
2912     for (i = 0; i < nelem/2; i++)
2913     {
2914         j = 0;
2915         while ((j < nr) &&
2916                gmx_strcasecmp(names[2*i], *(gnames[groups->grps[egcENER].nm_ind[j]])))
2917         {
2918             j++;
2919         }
2920         if (j == nr)
2921         {
2922             gmx_fatal(FARGS, "%s in %s is not an energy group\n",
2923                       names[2*i], option);
2924         }
2925         k = 0;
2926         while ((k < nr) &&
2927                gmx_strcasecmp(names[2*i+1], *(gnames[groups->grps[egcENER].nm_ind[k]])))
2928         {
2929             k++;
2930         }
2931         if (k == nr)
2932         {
2933             gmx_fatal(FARGS, "%s in %s is not an energy group\n",
2934                       names[2*i+1], option);
2935         }
2936         if ((j < nr) && (k < nr))
2937         {
2938             ir->opts.egp_flags[nr*j+k] |= flag;
2939             ir->opts.egp_flags[nr*k+j] |= flag;
2940             bSet = TRUE;
2941         }
2942     }
2943
2944     return bSet;
2945 }
2946
2947
2948 static void make_swap_groups(
2949         t_swapcoords *swap,
2950         char         *swapgname,
2951         char         *splitg0name,
2952         char         *splitg1name,
2953         char         *solgname,
2954         t_blocka     *grps,
2955         char        **gnames)
2956 {
2957     int   ig = -1, i = 0, j;
2958     char *splitg;
2959
2960
2961     /* Just a quick check here, more thorough checks are in mdrun */
2962     if (strcmp(splitg0name, splitg1name) == 0)
2963     {
2964         gmx_fatal(FARGS, "The split groups can not both be '%s'.", splitg0name);
2965     }
2966
2967     /* First get the swap group index atoms */
2968     ig        = search_string(swapgname, grps->nr, gnames);
2969     swap->nat = grps->index[ig+1] - grps->index[ig];
2970     if (swap->nat > 0)
2971     {
2972         fprintf(stderr, "Swap group '%s' contains %d atoms.\n", swapgname, swap->nat);
2973         snew(swap->ind, swap->nat);
2974         for (i = 0; i < swap->nat; i++)
2975         {
2976             swap->ind[i] = grps->a[grps->index[ig]+i];
2977         }
2978     }
2979     else
2980     {
2981         gmx_fatal(FARGS, "You defined an empty group of atoms for swapping.");
2982     }
2983
2984     /* Now do so for the split groups */
2985     for (j = 0; j < 2; j++)
2986     {
2987         if (j == 0)
2988         {
2989             splitg = splitg0name;
2990         }
2991         else
2992         {
2993             splitg = splitg1name;
2994         }
2995
2996         ig                 = search_string(splitg, grps->nr, gnames);
2997         swap->nat_split[j] = grps->index[ig+1] - grps->index[ig];
2998         if (swap->nat_split[j] > 0)
2999         {
3000             fprintf(stderr, "Split group %d '%s' contains %d atom%s.\n",
3001                     j, splitg, swap->nat_split[j], (swap->nat_split[j] > 1) ? "s" : "");
3002             snew(swap->ind_split[j], swap->nat_split[j]);
3003             for (i = 0; i < swap->nat_split[j]; i++)
3004             {
3005                 swap->ind_split[j][i] = grps->a[grps->index[ig]+i];
3006             }
3007         }
3008         else
3009         {
3010             gmx_fatal(FARGS, "Split group %d has to contain at least 1 atom!", j);
3011         }
3012     }
3013
3014     /* Now get the solvent group index atoms */
3015     ig            = search_string(solgname, grps->nr, gnames);
3016     swap->nat_sol = grps->index[ig+1] - grps->index[ig];
3017     if (swap->nat_sol > 0)
3018     {
3019         fprintf(stderr, "Solvent group '%s' contains %d atoms.\n", solgname, swap->nat_sol);
3020         snew(swap->ind_sol, swap->nat_sol);
3021         for (i = 0; i < swap->nat_sol; i++)
3022         {
3023             swap->ind_sol[i] = grps->a[grps->index[ig]+i];
3024         }
3025     }
3026     else
3027     {
3028         gmx_fatal(FARGS, "You defined an empty group of solvent. Cannot exchange ions.");
3029     }
3030 }
3031
3032
3033 void do_index(const char* mdparin, const char *ndx,
3034               gmx_mtop_t *mtop,
3035               gmx_bool bVerbose,
3036               t_inputrec *ir, rvec *v,
3037               warninp_t wi)
3038 {
3039     t_blocka     *grps;
3040     gmx_groups_t *groups;
3041     int           natoms;
3042     t_symtab     *symtab;
3043     t_atoms       atoms_all;
3044     char          warnbuf[STRLEN], **gnames;
3045     int           nr, ntcg, ntau_t, nref_t, nacc, nofg, nSA, nSA_points, nSA_time, nSA_temp;
3046     real          tau_min;
3047     int           nstcmin;
3048     int           nacg, nfreeze, nfrdim, nenergy, nvcm, nuser;
3049     char         *ptr1[MAXPTR], *ptr2[MAXPTR], *ptr3[MAXPTR];
3050     int           i, j, k, restnm;
3051     real          SAtime;
3052     gmx_bool      bExcl, bTable, bSetTCpar, bAnneal, bRest;
3053     int           nQMmethod, nQMbasis, nQMcharge, nQMmult, nbSH, nCASorb, nCASelec,
3054                   nSAon, nSAoff, nSAsteps, nQMg, nbOPT, nbTS;
3055     char          warn_buf[STRLEN];
3056
3057     if (bVerbose)
3058     {
3059         fprintf(stderr, "processing index file...\n");
3060     }
3061     debug_gmx();
3062     if (ndx == NULL)
3063     {
3064         snew(grps, 1);
3065         snew(grps->index, 1);
3066         snew(gnames, 1);
3067         atoms_all = gmx_mtop_global_atoms(mtop);
3068         analyse(&atoms_all, grps, &gnames, FALSE, TRUE);
3069         free_t_atoms(&atoms_all, FALSE);
3070     }
3071     else
3072     {
3073         grps = init_index(ndx, &gnames);
3074     }
3075
3076     groups = &mtop->groups;
3077     natoms = mtop->natoms;
3078     symtab = &mtop->symtab;
3079
3080     snew(groups->grpname, grps->nr+1);
3081
3082     for (i = 0; (i < grps->nr); i++)
3083     {
3084         groups->grpname[i] = put_symtab(symtab, gnames[i]);
3085     }
3086     groups->grpname[i] = put_symtab(symtab, "rest");
3087     restnm             = i;
3088     srenew(gnames, grps->nr+1);
3089     gnames[restnm]   = *(groups->grpname[i]);
3090     groups->ngrpname = grps->nr+1;
3091
3092     set_warning_line(wi, mdparin, -1);
3093
3094     ntau_t = str_nelem(is->tau_t, MAXPTR, ptr1);
3095     nref_t = str_nelem(is->ref_t, MAXPTR, ptr2);
3096     ntcg   = str_nelem(is->tcgrps, MAXPTR, ptr3);
3097     if ((ntau_t != ntcg) || (nref_t != ntcg))
3098     {
3099         gmx_fatal(FARGS, "Invalid T coupling input: %d groups, %d ref-t values and "
3100                   "%d tau-t values", ntcg, nref_t, ntau_t);
3101     }
3102
3103     bSetTCpar = (ir->etc || EI_SD(ir->eI) || ir->eI == eiBD || EI_TPI(ir->eI));
3104     do_numbering(natoms, groups, ntcg, ptr3, grps, gnames, egcTC,
3105                  restnm, bSetTCpar ? egrptpALL : egrptpALL_GENREST, bVerbose, wi);
3106     nr            = groups->grps[egcTC].nr;
3107     ir->opts.ngtc = nr;
3108     snew(ir->opts.nrdf, nr);
3109     snew(ir->opts.tau_t, nr);
3110     snew(ir->opts.ref_t, nr);
3111     if (ir->eI == eiBD && ir->bd_fric == 0)
3112     {
3113         fprintf(stderr, "bd-fric=0, so tau-t will be used as the inverse friction constant(s)\n");
3114     }
3115
3116     if (bSetTCpar)
3117     {
3118         if (nr != nref_t)
3119         {
3120             gmx_fatal(FARGS, "Not enough ref-t and tau-t values!");
3121         }
3122
3123         tau_min = 1e20;
3124         for (i = 0; (i < nr); i++)
3125         {
3126             ir->opts.tau_t[i] = strtod(ptr1[i], NULL);
3127             if ((ir->eI == eiBD || ir->eI == eiSD2) && ir->opts.tau_t[i] <= 0)
3128             {
3129                 sprintf(warn_buf, "With integrator %s tau-t should be larger than 0", ei_names[ir->eI]);
3130                 warning_error(wi, warn_buf);
3131             }
3132
3133             if (ir->etc != etcVRESCALE && ir->opts.tau_t[i] == 0)
3134             {
3135                 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.");
3136             }
3137
3138             if (ir->opts.tau_t[i] >= 0)
3139             {
3140                 tau_min = min(tau_min, ir->opts.tau_t[i]);
3141             }
3142         }
3143         if (ir->etc != etcNO && ir->nsttcouple == -1)
3144         {
3145             ir->nsttcouple = ir_optimal_nsttcouple(ir);
3146         }
3147
3148         if (EI_VV(ir->eI))
3149         {
3150             if ((ir->etc == etcNOSEHOOVER) && (ir->epc == epcBERENDSEN))
3151             {
3152                 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");
3153             }
3154             if ((ir->epc == epcMTTK) && (ir->etc > etcNO))
3155             {
3156                 if (ir->nstpcouple != ir->nsttcouple)
3157                 {
3158                     int mincouple = min(ir->nstpcouple, ir->nsttcouple);
3159                     ir->nstpcouple = ir->nsttcouple = mincouple;
3160                     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);
3161                     warning_note(wi, warn_buf);
3162                 }
3163             }
3164         }
3165         /* velocity verlet with averaged kinetic energy KE = 0.5*(v(t+1/2) - v(t-1/2)) is implemented
3166            primarily for testing purposes, and does not work with temperature coupling other than 1 */
3167
3168         if (ETC_ANDERSEN(ir->etc))
3169         {
3170             if (ir->nsttcouple != 1)
3171             {
3172                 ir->nsttcouple = 1;
3173                 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");
3174                 warning_note(wi, warn_buf);
3175             }
3176         }
3177         nstcmin = tcouple_min_integration_steps(ir->etc);
3178         if (nstcmin > 1)
3179         {
3180             if (tau_min/(ir->delta_t*ir->nsttcouple) < nstcmin)
3181             {
3182                 sprintf(warn_buf, "For proper integration of the %s thermostat, tau-t (%g) should be at least %d times larger than nsttcouple*dt (%g)",
3183                         ETCOUPLTYPE(ir->etc),
3184                         tau_min, nstcmin,
3185                         ir->nsttcouple*ir->delta_t);
3186                 warning(wi, warn_buf);
3187             }
3188         }
3189         for (i = 0; (i < nr); i++)
3190         {
3191             ir->opts.ref_t[i] = strtod(ptr2[i], NULL);
3192             if (ir->opts.ref_t[i] < 0)
3193             {
3194                 gmx_fatal(FARGS, "ref-t for group %d negative", i);
3195             }
3196         }
3197         /* set the lambda mc temperature to the md integrator temperature (which should be defined
3198            if we are in this conditional) if mc_temp is negative */
3199         if (ir->expandedvals->mc_temp < 0)
3200         {
3201             ir->expandedvals->mc_temp = ir->opts.ref_t[0]; /*for now, set to the first reft */
3202         }
3203     }
3204
3205     /* Simulated annealing for each group. There are nr groups */
3206     nSA = str_nelem(is->anneal, MAXPTR, ptr1);
3207     if (nSA == 1 && (ptr1[0][0] == 'n' || ptr1[0][0] == 'N'))
3208     {
3209         nSA = 0;
3210     }
3211     if (nSA > 0 && nSA != nr)
3212     {
3213         gmx_fatal(FARGS, "Not enough annealing values: %d (for %d groups)\n", nSA, nr);
3214     }
3215     else
3216     {
3217         snew(ir->opts.annealing, nr);
3218         snew(ir->opts.anneal_npoints, nr);
3219         snew(ir->opts.anneal_time, nr);
3220         snew(ir->opts.anneal_temp, nr);
3221         for (i = 0; i < nr; i++)
3222         {
3223             ir->opts.annealing[i]      = eannNO;
3224             ir->opts.anneal_npoints[i] = 0;
3225             ir->opts.anneal_time[i]    = NULL;
3226             ir->opts.anneal_temp[i]    = NULL;
3227         }
3228         if (nSA > 0)
3229         {
3230             bAnneal = FALSE;
3231             for (i = 0; i < nr; i++)
3232             {
3233                 if (ptr1[i][0] == 'n' || ptr1[i][0] == 'N')
3234                 {
3235                     ir->opts.annealing[i] = eannNO;
3236                 }
3237                 else if (ptr1[i][0] == 's' || ptr1[i][0] == 'S')
3238                 {
3239                     ir->opts.annealing[i] = eannSINGLE;
3240                     bAnneal               = TRUE;
3241                 }
3242                 else if (ptr1[i][0] == 'p' || ptr1[i][0] == 'P')
3243                 {
3244                     ir->opts.annealing[i] = eannPERIODIC;
3245                     bAnneal               = TRUE;
3246                 }
3247             }
3248             if (bAnneal)
3249             {
3250                 /* Read the other fields too */
3251                 nSA_points = str_nelem(is->anneal_npoints, MAXPTR, ptr1);
3252                 if (nSA_points != nSA)
3253                 {
3254                     gmx_fatal(FARGS, "Found %d annealing-npoints values for %d groups\n", nSA_points, nSA);
3255                 }
3256                 for (k = 0, i = 0; i < nr; i++)
3257                 {
3258                     ir->opts.anneal_npoints[i] = strtol(ptr1[i], NULL, 10);
3259                     if (ir->opts.anneal_npoints[i] == 1)
3260                     {
3261                         gmx_fatal(FARGS, "Please specify at least a start and an end point for annealing\n");
3262                     }
3263                     snew(ir->opts.anneal_time[i], ir->opts.anneal_npoints[i]);
3264                     snew(ir->opts.anneal_temp[i], ir->opts.anneal_npoints[i]);
3265                     k += ir->opts.anneal_npoints[i];
3266                 }
3267
3268                 nSA_time = str_nelem(is->anneal_time, MAXPTR, ptr1);
3269                 if (nSA_time != k)
3270                 {
3271                     gmx_fatal(FARGS, "Found %d annealing-time values, wanter %d\n", nSA_time, k);
3272                 }
3273                 nSA_temp = str_nelem(is->anneal_temp, MAXPTR, ptr2);
3274                 if (nSA_temp != k)
3275                 {
3276                     gmx_fatal(FARGS, "Found %d annealing-temp values, wanted %d\n", nSA_temp, k);
3277                 }
3278
3279                 for (i = 0, k = 0; i < nr; i++)
3280                 {
3281
3282                     for (j = 0; j < ir->opts.anneal_npoints[i]; j++)
3283                     {
3284                         ir->opts.anneal_time[i][j] = strtod(ptr1[k], NULL);
3285                         ir->opts.anneal_temp[i][j] = strtod(ptr2[k], NULL);
3286                         if (j == 0)
3287                         {
3288                             if (ir->opts.anneal_time[i][0] > (ir->init_t+GMX_REAL_EPS))
3289                             {
3290                                 gmx_fatal(FARGS, "First time point for annealing > init_t.\n");
3291                             }
3292                         }
3293                         else
3294                         {
3295                             /* j>0 */
3296                             if (ir->opts.anneal_time[i][j] < ir->opts.anneal_time[i][j-1])
3297                             {
3298                                 gmx_fatal(FARGS, "Annealing timepoints out of order: t=%f comes after t=%f\n",
3299                                           ir->opts.anneal_time[i][j], ir->opts.anneal_time[i][j-1]);
3300                             }
3301                         }
3302                         if (ir->opts.anneal_temp[i][j] < 0)
3303                         {
3304                             gmx_fatal(FARGS, "Found negative temperature in annealing: %f\n", ir->opts.anneal_temp[i][j]);
3305                         }
3306                         k++;
3307                     }
3308                 }
3309                 /* Print out some summary information, to make sure we got it right */
3310                 for (i = 0, k = 0; i < nr; i++)
3311                 {
3312                     if (ir->opts.annealing[i] != eannNO)
3313                     {
3314                         j = groups->grps[egcTC].nm_ind[i];
3315                         fprintf(stderr, "Simulated annealing for group %s: %s, %d timepoints\n",
3316                                 *(groups->grpname[j]), eann_names[ir->opts.annealing[i]],
3317                                 ir->opts.anneal_npoints[i]);
3318                         fprintf(stderr, "Time (ps)   Temperature (K)\n");
3319                         /* All terms except the last one */
3320                         for (j = 0; j < (ir->opts.anneal_npoints[i]-1); j++)
3321                         {
3322                             fprintf(stderr, "%9.1f      %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3323                         }
3324
3325                         /* Finally the last one */
3326                         j = ir->opts.anneal_npoints[i]-1;
3327                         if (ir->opts.annealing[i] == eannSINGLE)
3328                         {
3329                             fprintf(stderr, "%9.1f-     %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3330                         }
3331                         else
3332                         {
3333                             fprintf(stderr, "%9.1f      %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3334                             if (fabs(ir->opts.anneal_temp[i][j]-ir->opts.anneal_temp[i][0]) > GMX_REAL_EPS)
3335                             {
3336                                 warning_note(wi, "There is a temperature jump when your annealing loops back.\n");
3337                             }
3338                         }
3339                     }
3340                 }
3341             }
3342         }
3343     }
3344
3345     if (ir->ePull != epullNO)
3346     {
3347         make_pull_groups(ir->pull, is->pull_grp, grps, gnames);
3348
3349         make_pull_coords(ir->pull);
3350     }
3351
3352     if (ir->bRot)
3353     {
3354         make_rotation_groups(ir->rot, is->rot_grp, grps, gnames);
3355     }
3356
3357     if (ir->eSwapCoords != eswapNO)
3358     {
3359         make_swap_groups(ir->swap, swapgrp, splitgrp0, splitgrp1, solgrp, grps, gnames);
3360     }
3361
3362     nacc = str_nelem(is->acc, MAXPTR, ptr1);
3363     nacg = str_nelem(is->accgrps, MAXPTR, ptr2);
3364     if (nacg*DIM != nacc)
3365     {
3366         gmx_fatal(FARGS, "Invalid Acceleration input: %d groups and %d acc. values",
3367                   nacg, nacc);
3368     }
3369     do_numbering(natoms, groups, nacg, ptr2, grps, gnames, egcACC,
3370                  restnm, egrptpALL_GENREST, bVerbose, wi);
3371     nr = groups->grps[egcACC].nr;
3372     snew(ir->opts.acc, nr);
3373     ir->opts.ngacc = nr;
3374
3375     for (i = k = 0; (i < nacg); i++)
3376     {
3377         for (j = 0; (j < DIM); j++, k++)
3378         {
3379             ir->opts.acc[i][j] = strtod(ptr1[k], NULL);
3380         }
3381     }
3382     for (; (i < nr); i++)
3383     {
3384         for (j = 0; (j < DIM); j++)
3385         {
3386             ir->opts.acc[i][j] = 0;
3387         }
3388     }
3389
3390     nfrdim  = str_nelem(is->frdim, MAXPTR, ptr1);
3391     nfreeze = str_nelem(is->freeze, MAXPTR, ptr2);
3392     if (nfrdim != DIM*nfreeze)
3393     {
3394         gmx_fatal(FARGS, "Invalid Freezing input: %d groups and %d freeze values",
3395                   nfreeze, nfrdim);
3396     }
3397     do_numbering(natoms, groups, nfreeze, ptr2, grps, gnames, egcFREEZE,
3398                  restnm, egrptpALL_GENREST, bVerbose, wi);
3399     nr             = groups->grps[egcFREEZE].nr;
3400     ir->opts.ngfrz = nr;
3401     snew(ir->opts.nFreeze, nr);
3402     for (i = k = 0; (i < nfreeze); i++)
3403     {
3404         for (j = 0; (j < DIM); j++, k++)
3405         {
3406             ir->opts.nFreeze[i][j] = (gmx_strncasecmp(ptr1[k], "Y", 1) == 0);
3407             if (!ir->opts.nFreeze[i][j])
3408             {
3409                 if (gmx_strncasecmp(ptr1[k], "N", 1) != 0)
3410                 {
3411                     sprintf(warnbuf, "Please use Y(ES) or N(O) for freezedim only "
3412                             "(not %s)", ptr1[k]);
3413                     warning(wi, warn_buf);
3414                 }
3415             }
3416         }
3417     }
3418     for (; (i < nr); i++)
3419     {
3420         for (j = 0; (j < DIM); j++)
3421         {
3422             ir->opts.nFreeze[i][j] = 0;
3423         }
3424     }
3425
3426     nenergy = str_nelem(is->energy, MAXPTR, ptr1);
3427     do_numbering(natoms, groups, nenergy, ptr1, grps, gnames, egcENER,
3428                  restnm, egrptpALL_GENREST, bVerbose, wi);
3429     add_wall_energrps(groups, ir->nwall, symtab);
3430     ir->opts.ngener = groups->grps[egcENER].nr;
3431     nvcm            = str_nelem(is->vcm, MAXPTR, ptr1);
3432     bRest           =
3433         do_numbering(natoms, groups, nvcm, ptr1, grps, gnames, egcVCM,
3434                      restnm, nvcm == 0 ? egrptpALL_GENREST : egrptpPART, bVerbose, wi);
3435     if (bRest)
3436     {
3437         warning(wi, "Some atoms are not part of any center of mass motion removal group.\n"
3438                 "This may lead to artifacts.\n"
3439                 "In most cases one should use one group for the whole system.");
3440     }
3441
3442     /* Now we have filled the freeze struct, so we can calculate NRDF */
3443     calc_nrdf(mtop, ir, gnames);
3444
3445     if (v && NULL)
3446     {
3447         real fac, ntot = 0;
3448
3449         /* Must check per group! */
3450         for (i = 0; (i < ir->opts.ngtc); i++)
3451         {
3452             ntot += ir->opts.nrdf[i];
3453         }
3454         if (ntot != (DIM*natoms))
3455         {
3456             fac = sqrt(ntot/(DIM*natoms));
3457             if (bVerbose)
3458             {
3459                 fprintf(stderr, "Scaling velocities by a factor of %.3f to account for constraints\n"
3460                         "and removal of center of mass motion\n", fac);
3461             }
3462             for (i = 0; (i < natoms); i++)
3463             {
3464                 svmul(fac, v[i], v[i]);
3465             }
3466         }
3467     }
3468
3469     nuser = str_nelem(is->user1, MAXPTR, ptr1);
3470     do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcUser1,
3471                  restnm, egrptpALL_GENREST, bVerbose, wi);
3472     nuser = str_nelem(is->user2, MAXPTR, ptr1);
3473     do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcUser2,
3474                  restnm, egrptpALL_GENREST, bVerbose, wi);
3475     nuser = str_nelem(is->x_compressed_groups, MAXPTR, ptr1);
3476     do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcCompressedX,
3477                  restnm, egrptpONE, bVerbose, wi);
3478     nofg = str_nelem(is->orirefitgrp, MAXPTR, ptr1);
3479     do_numbering(natoms, groups, nofg, ptr1, grps, gnames, egcORFIT,
3480                  restnm, egrptpALL_GENREST, bVerbose, wi);
3481
3482     /* QMMM input processing */
3483     nQMg          = str_nelem(is->QMMM, MAXPTR, ptr1);
3484     nQMmethod     = str_nelem(is->QMmethod, MAXPTR, ptr2);
3485     nQMbasis      = str_nelem(is->QMbasis, MAXPTR, ptr3);
3486     if ((nQMmethod != nQMg) || (nQMbasis != nQMg))
3487     {
3488         gmx_fatal(FARGS, "Invalid QMMM input: %d groups %d basissets"
3489                   " and %d methods\n", nQMg, nQMbasis, nQMmethod);
3490     }
3491     /* group rest, if any, is always MM! */
3492     do_numbering(natoms, groups, nQMg, ptr1, grps, gnames, egcQMMM,
3493                  restnm, egrptpALL_GENREST, bVerbose, wi);
3494     nr            = nQMg; /*atoms->grps[egcQMMM].nr;*/
3495     ir->opts.ngQM = nQMg;
3496     snew(ir->opts.QMmethod, nr);
3497     snew(ir->opts.QMbasis, nr);
3498     for (i = 0; i < nr; i++)
3499     {
3500         /* input consists of strings: RHF CASSCF PM3 .. These need to be
3501          * converted to the corresponding enum in names.c
3502          */
3503         ir->opts.QMmethod[i] = search_QMstring(ptr2[i], eQMmethodNR,
3504                                                eQMmethod_names);
3505         ir->opts.QMbasis[i]  = search_QMstring(ptr3[i], eQMbasisNR,
3506                                                eQMbasis_names);
3507
3508     }
3509     nQMmult   = str_nelem(is->QMmult, MAXPTR, ptr1);
3510     nQMcharge = str_nelem(is->QMcharge, MAXPTR, ptr2);
3511     nbSH      = str_nelem(is->bSH, MAXPTR, ptr3);
3512     snew(ir->opts.QMmult, nr);
3513     snew(ir->opts.QMcharge, nr);
3514     snew(ir->opts.bSH, nr);
3515
3516     for (i = 0; i < nr; i++)
3517     {
3518         ir->opts.QMmult[i]   = strtol(ptr1[i], NULL, 10);
3519         ir->opts.QMcharge[i] = strtol(ptr2[i], NULL, 10);
3520         ir->opts.bSH[i]      = (gmx_strncasecmp(ptr3[i], "Y", 1) == 0);
3521     }
3522
3523     nCASelec  = str_nelem(is->CASelectrons, MAXPTR, ptr1);
3524     nCASorb   = str_nelem(is->CASorbitals, MAXPTR, ptr2);
3525     snew(ir->opts.CASelectrons, nr);
3526     snew(ir->opts.CASorbitals, nr);
3527     for (i = 0; i < nr; i++)
3528     {
3529         ir->opts.CASelectrons[i] = strtol(ptr1[i], NULL, 10);
3530         ir->opts.CASorbitals[i]  = strtol(ptr2[i], NULL, 10);
3531     }
3532     /* special optimization options */
3533
3534     nbOPT = str_nelem(is->bOPT, MAXPTR, ptr1);
3535     nbTS  = str_nelem(is->bTS, MAXPTR, ptr2);
3536     snew(ir->opts.bOPT, nr);
3537     snew(ir->opts.bTS, nr);
3538     for (i = 0; i < nr; i++)
3539     {
3540         ir->opts.bOPT[i] = (gmx_strncasecmp(ptr1[i], "Y", 1) == 0);
3541         ir->opts.bTS[i]  = (gmx_strncasecmp(ptr2[i], "Y", 1) == 0);
3542     }
3543     nSAon     = str_nelem(is->SAon, MAXPTR, ptr1);
3544     nSAoff    = str_nelem(is->SAoff, MAXPTR, ptr2);
3545     nSAsteps  = str_nelem(is->SAsteps, MAXPTR, ptr3);
3546     snew(ir->opts.SAon, nr);
3547     snew(ir->opts.SAoff, nr);
3548     snew(ir->opts.SAsteps, nr);
3549
3550     for (i = 0; i < nr; i++)
3551     {
3552         ir->opts.SAon[i]    = strtod(ptr1[i], NULL);
3553         ir->opts.SAoff[i]   = strtod(ptr2[i], NULL);
3554         ir->opts.SAsteps[i] = strtol(ptr3[i], NULL, 10);
3555     }
3556     /* end of QMMM input */
3557
3558     if (bVerbose)
3559     {
3560         for (i = 0; (i < egcNR); i++)
3561         {
3562             fprintf(stderr, "%-16s has %d element(s):", gtypes[i], groups->grps[i].nr);
3563             for (j = 0; (j < groups->grps[i].nr); j++)
3564             {
3565                 fprintf(stderr, " %s", *(groups->grpname[groups->grps[i].nm_ind[j]]));
3566             }
3567             fprintf(stderr, "\n");
3568         }
3569     }
3570
3571     nr = groups->grps[egcENER].nr;
3572     snew(ir->opts.egp_flags, nr*nr);
3573
3574     bExcl = do_egp_flag(ir, groups, "energygrp-excl", is->egpexcl, EGP_EXCL);
3575     if (bExcl && ir->cutoff_scheme == ecutsVERLET)
3576     {
3577         warning_error(wi, "Energy group exclusions are not (yet) implemented for the Verlet scheme");
3578     }
3579     if (bExcl && EEL_FULL(ir->coulombtype))
3580     {
3581         warning(wi, "Can not exclude the lattice Coulomb energy between energy groups");
3582     }
3583
3584     bTable = do_egp_flag(ir, groups, "energygrp-table", is->egptable, EGP_TABLE);
3585     if (bTable && !(ir->vdwtype == evdwUSER) &&
3586         !(ir->coulombtype == eelUSER) && !(ir->coulombtype == eelPMEUSER) &&
3587         !(ir->coulombtype == eelPMEUSERSWITCH))
3588     {
3589         gmx_fatal(FARGS, "Can only have energy group pair tables in combination with user tables for VdW and/or Coulomb");
3590     }
3591
3592     decode_cos(is->efield_x, &(ir->ex[XX]));
3593     decode_cos(is->efield_xt, &(ir->et[XX]));
3594     decode_cos(is->efield_y, &(ir->ex[YY]));
3595     decode_cos(is->efield_yt, &(ir->et[YY]));
3596     decode_cos(is->efield_z, &(ir->ex[ZZ]));
3597     decode_cos(is->efield_zt, &(ir->et[ZZ]));
3598
3599     if (ir->bAdress)
3600     {
3601         do_adress_index(ir->adress, groups, gnames, &(ir->opts), wi);
3602     }
3603
3604     for (i = 0; (i < grps->nr); i++)
3605     {
3606         sfree(gnames[i]);
3607     }
3608     sfree(gnames);
3609     done_blocka(grps);
3610     sfree(grps);
3611
3612 }
3613
3614
3615
3616 static void check_disre(gmx_mtop_t *mtop)
3617 {
3618     gmx_ffparams_t *ffparams;
3619     t_functype     *functype;
3620     t_iparams      *ip;
3621     int             i, ndouble, ftype;
3622     int             label, old_label;
3623
3624     if (gmx_mtop_ftype_count(mtop, F_DISRES) > 0)
3625     {
3626         ffparams  = &mtop->ffparams;
3627         functype  = ffparams->functype;
3628         ip        = ffparams->iparams;
3629         ndouble   = 0;
3630         old_label = -1;
3631         for (i = 0; i < ffparams->ntypes; i++)
3632         {
3633             ftype = functype[i];
3634             if (ftype == F_DISRES)
3635             {
3636                 label = ip[i].disres.label;
3637                 if (label == old_label)
3638                 {
3639                     fprintf(stderr, "Distance restraint index %d occurs twice\n", label);
3640                     ndouble++;
3641                 }
3642                 old_label = label;
3643             }
3644         }
3645         if (ndouble > 0)
3646         {
3647             gmx_fatal(FARGS, "Found %d double distance restraint indices,\n"
3648                       "probably the parameters for multiple pairs in one restraint "
3649                       "are not identical\n", ndouble);
3650         }
3651     }
3652 }
3653
3654 static gmx_bool absolute_reference(t_inputrec *ir, gmx_mtop_t *sys,
3655                                    gmx_bool posres_only,
3656                                    ivec AbsRef)
3657 {
3658     int                  d, g, i;
3659     gmx_mtop_ilistloop_t iloop;
3660     t_ilist             *ilist;
3661     int                  nmol;
3662     t_iparams           *pr;
3663
3664     clear_ivec(AbsRef);
3665
3666     if (!posres_only)
3667     {
3668         /* Check the COM */
3669         for (d = 0; d < DIM; d++)
3670         {
3671             AbsRef[d] = (d < ndof_com(ir) ? 0 : 1);
3672         }
3673         /* Check for freeze groups */
3674         for (g = 0; g < ir->opts.ngfrz; g++)
3675         {
3676             for (d = 0; d < DIM; d++)
3677             {
3678                 if (ir->opts.nFreeze[g][d] != 0)
3679                 {
3680                     AbsRef[d] = 1;
3681                 }
3682             }
3683         }
3684     }
3685
3686     /* Check for position restraints */
3687     iloop = gmx_mtop_ilistloop_init(sys);
3688     while (gmx_mtop_ilistloop_next(iloop, &ilist, &nmol))
3689     {
3690         if (nmol > 0 &&
3691             (AbsRef[XX] == 0 || AbsRef[YY] == 0 || AbsRef[ZZ] == 0))
3692         {
3693             for (i = 0; i < ilist[F_POSRES].nr; i += 2)
3694             {
3695                 pr = &sys->ffparams.iparams[ilist[F_POSRES].iatoms[i]];
3696                 for (d = 0; d < DIM; d++)
3697                 {
3698                     if (pr->posres.fcA[d] != 0)
3699                     {
3700                         AbsRef[d] = 1;
3701                     }
3702                 }
3703             }
3704             for (i = 0; i < ilist[F_FBPOSRES].nr; i += 2)
3705             {
3706                 /* Check for flat-bottom posres */
3707                 pr = &sys->ffparams.iparams[ilist[F_FBPOSRES].iatoms[i]];
3708                 if (pr->fbposres.k != 0)
3709                 {
3710                     switch (pr->fbposres.geom)
3711                     {
3712                         case efbposresSPHERE:
3713                             AbsRef[XX] = AbsRef[YY] = AbsRef[ZZ] = 1;
3714                             break;
3715                         case efbposresCYLINDER:
3716                             AbsRef[XX] = AbsRef[YY] = 1;
3717                             break;
3718                         case efbposresX: /* d=XX */
3719                         case efbposresY: /* d=YY */
3720                         case efbposresZ: /* d=ZZ */
3721                             d         = pr->fbposres.geom - efbposresX;
3722                             AbsRef[d] = 1;
3723                             break;
3724                         default:
3725                             gmx_fatal(FARGS, " Invalid geometry for flat-bottom position restraint.\n"
3726                                       "Expected nr between 1 and %d. Found %d\n", efbposresNR-1,
3727                                       pr->fbposres.geom);
3728                     }
3729                 }
3730             }
3731         }
3732     }
3733
3734     return (AbsRef[XX] != 0 && AbsRef[YY] != 0 && AbsRef[ZZ] != 0);
3735 }
3736
3737 static void
3738 check_combination_rule_differences(const gmx_mtop_t *mtop, int state,
3739                                    gmx_bool *bC6ParametersWorkWithGeometricRules,
3740                                    gmx_bool *bC6ParametersWorkWithLBRules,
3741                                    gmx_bool *bLBRulesPossible)
3742 {
3743     int           ntypes, tpi, tpj, thisLBdiff, thisgeomdiff;
3744     int          *typecount;
3745     real          tol;
3746     double        geometricdiff, LBdiff;
3747     double        c6i, c6j, c12i, c12j;
3748     double        c6, c6_geometric, c6_LB;
3749     double        sigmai, sigmaj, epsi, epsj;
3750     gmx_bool      bCanDoLBRules, bCanDoGeometricRules;
3751     const char   *ptr;
3752
3753     /* A tolerance of 1e-5 seems reasonable for (possibly hand-typed)
3754      * force-field floating point parameters.
3755      */
3756     tol = 1e-5;
3757     ptr = getenv("GMX_LJCOMB_TOL");
3758     if (ptr != NULL)
3759     {
3760         double dbl;
3761
3762         sscanf(ptr, "%lf", &dbl);
3763         tol = dbl;
3764     }
3765
3766     *bC6ParametersWorkWithLBRules         = TRUE;
3767     *bC6ParametersWorkWithGeometricRules  = TRUE;
3768     bCanDoLBRules                         = TRUE;
3769     bCanDoGeometricRules                  = TRUE;
3770     ntypes                                = mtop->ffparams.atnr;
3771     snew(typecount, ntypes);
3772     gmx_mtop_count_atomtypes(mtop, state, typecount);
3773     geometricdiff           = LBdiff = 0.0;
3774     *bLBRulesPossible       = TRUE;
3775     for (tpi = 0; tpi < ntypes; ++tpi)
3776     {
3777         c6i  = mtop->ffparams.iparams[(ntypes + 1) * tpi].lj.c6;
3778         c12i = mtop->ffparams.iparams[(ntypes + 1) * tpi].lj.c12;
3779         for (tpj = tpi; tpj < ntypes; ++tpj)
3780         {
3781             c6j          = mtop->ffparams.iparams[(ntypes + 1) * tpj].lj.c6;
3782             c12j         = mtop->ffparams.iparams[(ntypes + 1) * tpj].lj.c12;
3783             c6           = mtop->ffparams.iparams[ntypes * tpi + tpj].lj.c6;
3784             c6_geometric = sqrt(c6i * c6j);
3785             if (!gmx_numzero(c6_geometric))
3786             {
3787                 if (!gmx_numzero(c12i) && !gmx_numzero(c12j))
3788                 {
3789                     sigmai   = pow(c12i / c6i, 1.0/6.0);
3790                     sigmaj   = pow(c12j / c6j, 1.0/6.0);
3791                     epsi     = c6i * c6i /(4.0 * c12i);
3792                     epsj     = c6j * c6j /(4.0 * c12j);
3793                     c6_LB    = 4.0 * pow(epsi * epsj, 1.0/2.0) * pow(0.5 * (sigmai + sigmaj), 6);
3794                 }
3795                 else
3796                 {
3797                     *bLBRulesPossible = FALSE;
3798                     c6_LB             = c6_geometric;
3799                 }
3800                 bCanDoLBRules = gmx_within_tol(c6_LB, c6, tol);
3801             }
3802
3803             if (FALSE == bCanDoLBRules)
3804             {
3805                 *bC6ParametersWorkWithLBRules = FALSE;
3806             }
3807
3808             bCanDoGeometricRules = gmx_within_tol(c6_geometric, c6, tol);
3809
3810             if (FALSE == bCanDoGeometricRules)
3811             {
3812                 *bC6ParametersWorkWithGeometricRules = FALSE;
3813             }
3814         }
3815     }
3816     sfree(typecount);
3817 }
3818
3819 static void
3820 check_combination_rules(const t_inputrec *ir, const gmx_mtop_t *mtop,
3821                         warninp_t wi)
3822 {
3823     char     err_buf[256];
3824     gmx_bool bLBRulesPossible, bC6ParametersWorkWithGeometricRules, bC6ParametersWorkWithLBRules;
3825
3826     check_combination_rule_differences(mtop, 0,
3827                                        &bC6ParametersWorkWithGeometricRules,
3828                                        &bC6ParametersWorkWithLBRules,
3829                                        &bLBRulesPossible);
3830     if (ir->ljpme_combination_rule == eljpmeLB)
3831     {
3832         if (FALSE == bC6ParametersWorkWithLBRules || FALSE == bLBRulesPossible)
3833         {
3834             warning(wi, "You are using arithmetic-geometric combination rules "
3835                     "in LJ-PME, but your non-bonded C6 parameters do not "
3836                     "follow these rules.");
3837         }
3838     }
3839     else
3840     {
3841         if (FALSE == bC6ParametersWorkWithGeometricRules)
3842         {
3843             if (ir->eDispCorr != edispcNO)
3844             {
3845                 warning_note(wi, "You are using geometric combination rules in "
3846                              "LJ-PME, but your non-bonded C6 parameters do "
3847                              "not follow these rules. "
3848                              "If your force field uses Lorentz-Berthelot combination rules this "
3849                              "will introduce small errors in the forces and energies in "
3850                              "your simulations. Dispersion correction will correct total "
3851                              "energy and/or pressure, but not forces or surface tensions. "
3852                              "Please check the LJ-PME section in the manual "
3853                              "before proceeding further.");
3854             }
3855             else
3856             {
3857                 warning_note(wi, "You are using geometric combination rules in "
3858                              "LJ-PME, but your non-bonded C6 parameters do "
3859                              "not follow these rules. "
3860                              "If your force field uses Lorentz-Berthelot combination rules this "
3861                              "will introduce small errors in the forces and energies in "
3862                              "your simulations. Consider using dispersion correction "
3863                              "for the total energy and pressure. "
3864                              "Please check the LJ-PME section in the manual "
3865                              "before proceeding further.");
3866             }
3867         }
3868     }
3869 }
3870
3871 void triple_check(const char *mdparin, t_inputrec *ir, gmx_mtop_t *sys,
3872                   warninp_t wi)
3873 {
3874     char                      err_buf[256];
3875     int                       i, m, c, nmol, npct;
3876     gmx_bool                  bCharge, bAcc;
3877     real                      gdt_max, *mgrp, mt;
3878     rvec                      acc;
3879     gmx_mtop_atomloop_block_t aloopb;
3880     gmx_mtop_atomloop_all_t   aloop;
3881     t_atom                   *atom;
3882     ivec                      AbsRef;
3883     char                      warn_buf[STRLEN];
3884
3885     set_warning_line(wi, mdparin, -1);
3886
3887     if (EI_DYNAMICS(ir->eI) && !EI_SD(ir->eI) && ir->eI != eiBD &&
3888         ir->comm_mode == ecmNO &&
3889         !(absolute_reference(ir, sys, FALSE, AbsRef) || ir->nsteps <= 10))
3890     {
3891         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");
3892     }
3893
3894     /* Check for pressure coupling with absolute position restraints */
3895     if (ir->epc != epcNO && ir->refcoord_scaling == erscNO)
3896     {
3897         absolute_reference(ir, sys, TRUE, AbsRef);
3898         {
3899             for (m = 0; m < DIM; m++)
3900             {
3901                 if (AbsRef[m] && norm2(ir->compress[m]) > 0)
3902                 {
3903                     warning(wi, "You are using pressure coupling with absolute position restraints, this will give artifacts. Use the refcoord_scaling option.");
3904                     break;
3905                 }
3906             }
3907         }
3908     }
3909
3910     bCharge = FALSE;
3911     aloopb  = gmx_mtop_atomloop_block_init(sys);
3912     while (gmx_mtop_atomloop_block_next(aloopb, &atom, &nmol))
3913     {
3914         if (atom->q != 0 || atom->qB != 0)
3915         {
3916             bCharge = TRUE;
3917         }
3918     }
3919
3920     if (!bCharge)
3921     {
3922         if (EEL_FULL(ir->coulombtype))
3923         {
3924             sprintf(err_buf,
3925                     "You are using full electrostatics treatment %s for a system without charges.\n"
3926                     "This costs a lot of performance for just processing zeros, consider using %s instead.\n",
3927                     EELTYPE(ir->coulombtype), EELTYPE(eelCUT));
3928             warning(wi, err_buf);
3929         }
3930     }
3931     else
3932     {
3933         if (ir->coulombtype == eelCUT && ir->rcoulomb > 0 && !ir->implicit_solvent)
3934         {
3935             sprintf(err_buf,
3936                     "You are using a plain Coulomb cut-off, which might produce artifacts.\n"
3937                     "You might want to consider using %s electrostatics.\n",
3938                     EELTYPE(eelPME));
3939             warning_note(wi, err_buf);
3940         }
3941     }
3942
3943     /* Check if combination rules used in LJ-PME are the same as in the force field */
3944     if (EVDW_PME(ir->vdwtype))
3945     {
3946         check_combination_rules(ir, sys, wi);
3947     }
3948
3949     /* Generalized reaction field */
3950     if (ir->opts.ngtc == 0)
3951     {
3952         sprintf(err_buf, "No temperature coupling while using coulombtype %s",
3953                 eel_names[eelGRF]);
3954         CHECK(ir->coulombtype == eelGRF);
3955     }
3956     else
3957     {
3958         sprintf(err_buf, "When using coulombtype = %s"
3959                 " ref-t for temperature coupling should be > 0",
3960                 eel_names[eelGRF]);
3961         CHECK((ir->coulombtype == eelGRF) && (ir->opts.ref_t[0] <= 0));
3962     }
3963
3964     if (ir->eI == eiSD1 &&
3965         (gmx_mtop_ftype_count(sys, F_CONSTR) > 0 ||
3966          gmx_mtop_ftype_count(sys, F_SETTLE) > 0))
3967     {
3968         sprintf(warn_buf, "With constraints integrator %s is less accurate, consider using %s instead", ei_names[ir->eI], ei_names[eiSD2]);
3969         warning_note(wi, warn_buf);
3970     }
3971
3972     bAcc = FALSE;
3973     for (i = 0; (i < sys->groups.grps[egcACC].nr); i++)
3974     {
3975         for (m = 0; (m < DIM); m++)
3976         {
3977             if (fabs(ir->opts.acc[i][m]) > 1e-6)
3978             {
3979                 bAcc = TRUE;
3980             }
3981         }
3982     }
3983     if (bAcc)
3984     {
3985         clear_rvec(acc);
3986         snew(mgrp, sys->groups.grps[egcACC].nr);
3987         aloop = gmx_mtop_atomloop_all_init(sys);
3988         while (gmx_mtop_atomloop_all_next(aloop, &i, &atom))
3989         {
3990             mgrp[ggrpnr(&sys->groups, egcACC, i)] += atom->m;
3991         }
3992         mt = 0.0;
3993         for (i = 0; (i < sys->groups.grps[egcACC].nr); i++)
3994         {
3995             for (m = 0; (m < DIM); m++)
3996             {
3997                 acc[m] += ir->opts.acc[i][m]*mgrp[i];
3998             }
3999             mt += mgrp[i];
4000         }
4001         for (m = 0; (m < DIM); m++)
4002         {
4003             if (fabs(acc[m]) > 1e-6)
4004             {
4005                 const char *dim[DIM] = { "X", "Y", "Z" };
4006                 fprintf(stderr,
4007                         "Net Acceleration in %s direction, will %s be corrected\n",
4008                         dim[m], ir->nstcomm != 0 ? "" : "not");
4009                 if (ir->nstcomm != 0 && m < ndof_com(ir))
4010                 {
4011                     acc[m] /= mt;
4012                     for (i = 0; (i < sys->groups.grps[egcACC].nr); i++)
4013                     {
4014                         ir->opts.acc[i][m] -= acc[m];
4015                     }
4016                 }
4017             }
4018         }
4019         sfree(mgrp);
4020     }
4021
4022     if (ir->efep != efepNO && ir->fepvals->sc_alpha != 0 &&
4023         !gmx_within_tol(sys->ffparams.reppow, 12.0, 10*GMX_DOUBLE_EPS))
4024     {
4025         gmx_fatal(FARGS, "Soft-core interactions are only supported with VdW repulsion power 12");
4026     }
4027
4028     if (ir->ePull != epullNO)
4029     {
4030         gmx_bool bPullAbsoluteRef;
4031
4032         bPullAbsoluteRef = FALSE;
4033         for (i = 0; i < ir->pull->ncoord; i++)
4034         {
4035             bPullAbsoluteRef = bPullAbsoluteRef ||
4036                 ir->pull->coord[i].group[0] == 0 ||
4037                 ir->pull->coord[i].group[1] == 0;
4038         }
4039         if (bPullAbsoluteRef)
4040         {
4041             absolute_reference(ir, sys, FALSE, AbsRef);
4042             for (m = 0; m < DIM; m++)
4043             {
4044                 if (ir->pull->dim[m] && !AbsRef[m])
4045                 {
4046                     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.");
4047                     break;
4048                 }
4049             }
4050         }
4051
4052         if (ir->pull->eGeom == epullgDIRPBC)
4053         {
4054             for (i = 0; i < 3; i++)
4055             {
4056                 for (m = 0; m <= i; m++)
4057                 {
4058                     if ((ir->epc != epcNO && ir->compress[i][m] != 0) ||
4059                         ir->deform[i][m] != 0)
4060                     {
4061                         for (c = 0; c < ir->pull->ncoord; c++)
4062                         {
4063                             if (ir->pull->coord[c].vec[m] != 0)
4064                             {
4065                                 gmx_fatal(FARGS, "Can not have dynamic box while using pull geometry '%s' (dim %c)", EPULLGEOM(ir->pull->eGeom), 'x'+m);
4066                             }
4067                         }
4068                     }
4069                 }
4070             }
4071         }
4072     }
4073
4074     check_disre(sys);
4075 }
4076
4077 void double_check(t_inputrec *ir, matrix box, gmx_bool bConstr, warninp_t wi)
4078 {
4079     real        min_size;
4080     gmx_bool    bTWIN;
4081     char        warn_buf[STRLEN];
4082     const char *ptr;
4083
4084     ptr = check_box(ir->ePBC, box);
4085     if (ptr)
4086     {
4087         warning_error(wi, ptr);
4088     }
4089
4090     if (bConstr && ir->eConstrAlg == econtSHAKE)
4091     {
4092         if (ir->shake_tol <= 0.0)
4093         {
4094             sprintf(warn_buf, "ERROR: shake-tol must be > 0 instead of %g\n",
4095                     ir->shake_tol);
4096             warning_error(wi, warn_buf);
4097         }
4098
4099         if (IR_TWINRANGE(*ir) && ir->nstlist > 1)
4100         {
4101             sprintf(warn_buf, "With twin-range cut-off's and SHAKE the virial and the pressure are incorrect.");
4102             if (ir->epc == epcNO)
4103             {
4104                 warning(wi, warn_buf);
4105             }
4106             else
4107             {
4108                 warning_error(wi, warn_buf);
4109             }
4110         }
4111     }
4112
4113     if ( (ir->eConstrAlg == econtLINCS) && bConstr)
4114     {
4115         /* If we have Lincs constraints: */
4116         if (ir->eI == eiMD && ir->etc == etcNO &&
4117             ir->eConstrAlg == econtLINCS && ir->nLincsIter == 1)
4118         {
4119             sprintf(warn_buf, "For energy conservation with LINCS, lincs_iter should be 2 or larger.\n");
4120             warning_note(wi, warn_buf);
4121         }
4122
4123         if ((ir->eI == eiCG || ir->eI == eiLBFGS) && (ir->nProjOrder < 8))
4124         {
4125             sprintf(warn_buf, "For accurate %s with LINCS constraints, lincs-order should be 8 or more.", ei_names[ir->eI]);
4126             warning_note(wi, warn_buf);
4127         }
4128         if (ir->epc == epcMTTK)
4129         {
4130             warning_error(wi, "MTTK not compatible with lincs -- use shake instead.");
4131         }
4132     }
4133
4134     if (ir->LincsWarnAngle > 90.0)
4135     {
4136         sprintf(warn_buf, "lincs-warnangle can not be larger than 90 degrees, setting it to 90.\n");
4137         warning(wi, warn_buf);
4138         ir->LincsWarnAngle = 90.0;
4139     }
4140
4141     if (ir->ePBC != epbcNONE)
4142     {
4143         if (ir->nstlist == 0)
4144         {
4145             warning(wi, "With nstlist=0 atoms are only put into the box at step 0, therefore drifting atoms might cause the simulation to crash.");
4146         }
4147         bTWIN = (ir->rlistlong > ir->rlist);
4148         if (ir->ns_type == ensGRID)
4149         {
4150             if (sqr(ir->rlistlong) >= max_cutoff2(ir->ePBC, box))
4151             {
4152                 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",
4153                         bTWIN ? (ir->rcoulomb == ir->rlistlong ? "rcoulomb" : "rvdw") : "rlist");
4154                 warning_error(wi, warn_buf);
4155             }
4156         }
4157         else
4158         {
4159             min_size = min(box[XX][XX], min(box[YY][YY], box[ZZ][ZZ]));
4160             if (2*ir->rlistlong >= min_size)
4161             {
4162                 sprintf(warn_buf, "ERROR: One of the box lengths is smaller than twice the cut-off length. Increase the box size or decrease rlist.");
4163                 warning_error(wi, warn_buf);
4164                 if (TRICLINIC(box))
4165                 {
4166                     fprintf(stderr, "Grid search might allow larger cut-off's than simple search with triclinic boxes.");
4167                 }
4168             }
4169         }
4170     }
4171 }
4172
4173 void check_chargegroup_radii(const gmx_mtop_t *mtop, const t_inputrec *ir,
4174                              rvec *x,
4175                              warninp_t wi)
4176 {
4177     real rvdw1, rvdw2, rcoul1, rcoul2;
4178     char warn_buf[STRLEN];
4179
4180     calc_chargegroup_radii(mtop, x, &rvdw1, &rvdw2, &rcoul1, &rcoul2);
4181
4182     if (rvdw1 > 0)
4183     {
4184         printf("Largest charge group radii for Van der Waals: %5.3f, %5.3f nm\n",
4185                rvdw1, rvdw2);
4186     }
4187     if (rcoul1 > 0)
4188     {
4189         printf("Largest charge group radii for Coulomb:       %5.3f, %5.3f nm\n",
4190                rcoul1, rcoul2);
4191     }
4192
4193     if (ir->rlist > 0)
4194     {
4195         if (rvdw1  + rvdw2  > ir->rlist ||
4196             rcoul1 + rcoul2 > ir->rlist)
4197         {
4198             sprintf(warn_buf,
4199                     "The sum of the two largest charge group radii (%f) "
4200                     "is larger than rlist (%f)\n",
4201                     max(rvdw1+rvdw2, rcoul1+rcoul2), ir->rlist);
4202             warning(wi, warn_buf);
4203         }
4204         else
4205         {
4206             /* Here we do not use the zero at cut-off macro,
4207              * since user defined interactions might purposely
4208              * not be zero at the cut-off.
4209              */
4210             if ((EVDW_IS_ZERO_AT_CUTOFF(ir->vdwtype) ||
4211                  ir->vdw_modifier != eintmodNONE) &&
4212                 rvdw1 + rvdw2 > ir->rlistlong - ir->rvdw)
4213             {
4214                 sprintf(warn_buf, "The sum of the two largest charge group "
4215                         "radii (%f) is larger than %s (%f) - rvdw (%f).\n"
4216                         "With exact cut-offs, better performance can be "
4217                         "obtained with cutoff-scheme = %s, because it "
4218                         "does not use charge groups at all.",
4219                         rvdw1+rvdw2,
4220                         ir->rlistlong > ir->rlist ? "rlistlong" : "rlist",
4221                         ir->rlistlong, ir->rvdw,
4222                         ecutscheme_names[ecutsVERLET]);
4223                 if (ir_NVE(ir))
4224                 {
4225                     warning(wi, warn_buf);
4226                 }
4227                 else
4228                 {
4229                     warning_note(wi, warn_buf);
4230                 }
4231             }
4232             if ((EEL_IS_ZERO_AT_CUTOFF(ir->coulombtype) ||
4233                  ir->coulomb_modifier != eintmodNONE) &&
4234                 rcoul1 + rcoul2 > ir->rlistlong - ir->rcoulomb)
4235             {
4236                 sprintf(warn_buf, "The sum of the two largest charge group radii (%f) is larger than %s (%f) - rcoulomb (%f).\n"
4237                         "With exact cut-offs, better performance can be obtained with cutoff-scheme = %s, because it does not use charge groups at all.",
4238                         rcoul1+rcoul2,
4239                         ir->rlistlong > ir->rlist ? "rlistlong" : "rlist",
4240                         ir->rlistlong, ir->rcoulomb,
4241                         ecutscheme_names[ecutsVERLET]);
4242                 if (ir_NVE(ir))
4243                 {
4244                     warning(wi, warn_buf);
4245                 }
4246                 else
4247                 {
4248                     warning_note(wi, warn_buf);
4249                 }
4250             }
4251         }
4252     }
4253 }