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