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