src/gromacs/gmxpreprocess/readir.c

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