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