src/gromacs/gmxpreprocess/readir.c

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