src/gromacs/gmxpreprocess/readir.cpp

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