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