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     CTYPE ("Periodic boundary conditions: xyz, no, xy");
1870     EETYPE("pbc",         ir->ePBC,       epbc_names);
1871     EETYPE("periodic-molecules", ir->bPeriodicMols, yesno_names);
1872     CTYPE ("Allowed energy error due to the Verlet buffer in kJ/mol/ps per atom,");
1873     CTYPE ("a value of -1 means: use rlist");
1874     RTYPE("verlet-buffer-tolerance", ir->verletbuf_tol,    0.005);
1875     CTYPE ("nblist cut-off");
1876     RTYPE ("rlist",   ir->rlist,  1.0);
1877     CTYPE ("long-range cut-off for switched potentials");
1878     RTYPE ("rlistlong",   ir->rlistlong,  -1);
1879     ITYPE ("nstcalclr",   ir->nstcalclr,  -1);
1880
1881     /* Electrostatics */
1882     CCTYPE ("OPTIONS FOR ELECTROSTATICS AND VDW");
1883     CTYPE ("Method for doing electrostatics");
1884     EETYPE("coulombtype", ir->coulombtype,    eel_names);
1885     EETYPE("coulomb-modifier",    ir->coulomb_modifier,    eintmod_names);
1886     CTYPE ("cut-off lengths");
1887     RTYPE ("rcoulomb-switch", ir->rcoulomb_switch,    0.0);
1888     RTYPE ("rcoulomb",    ir->rcoulomb,   1.0);
1889     CTYPE ("Relative dielectric constant for the medium and the reaction field");
1890     RTYPE ("epsilon-r",   ir->epsilon_r,  1.0);
1891     RTYPE ("epsilon-rf",  ir->epsilon_rf, 0.0);
1892     CTYPE ("Method for doing Van der Waals");
1893     EETYPE("vdw-type",    ir->vdwtype,    evdw_names);
1894     EETYPE("vdw-modifier",    ir->vdw_modifier,    eintmod_names);
1895     CTYPE ("cut-off lengths");
1896     RTYPE ("rvdw-switch", ir->rvdw_switch,    0.0);
1897     RTYPE ("rvdw",    ir->rvdw,   1.0);
1898     CTYPE ("Apply long range dispersion corrections for Energy and Pressure");
1899     EETYPE("DispCorr",    ir->eDispCorr,  edispc_names);
1900     CTYPE ("Extension of the potential lookup tables beyond the cut-off");
1901     RTYPE ("table-extension", ir->tabext, 1.0);
1902     CTYPE ("Separate tables between energy group pairs");
1903     STYPE ("energygrp-table", is->egptable,   NULL);
1904     CTYPE ("Spacing for the PME/PPPM FFT grid");
1905     RTYPE ("fourierspacing", ir->fourier_spacing, 0.12);
1906     CTYPE ("FFT grid size, when a value is 0 fourierspacing will be used");
1907     ITYPE ("fourier-nx",  ir->nkx,         0);
1908     ITYPE ("fourier-ny",  ir->nky,         0);
1909     ITYPE ("fourier-nz",  ir->nkz,         0);
1910     CTYPE ("EWALD/PME/PPPM parameters");
1911     ITYPE ("pme-order",   ir->pme_order,   4);
1912     RTYPE ("ewald-rtol",  ir->ewald_rtol, 0.00001);
1913     RTYPE ("ewald-rtol-lj", ir->ewald_rtol_lj, 0.001);
1914     EETYPE("lj-pme-comb-rule", ir->ljpme_combination_rule, eljpme_names);
1915     EETYPE("ewald-geometry", ir->ewald_geometry, eewg_names);
1916     RTYPE ("epsilon-surface", ir->epsilon_surface, 0.0);
1917
1918     CCTYPE("IMPLICIT SOLVENT ALGORITHM");
1919     EETYPE("implicit-solvent", ir->implicit_solvent, eis_names);
1920
1921     CCTYPE ("GENERALIZED BORN ELECTROSTATICS");
1922     CTYPE ("Algorithm for calculating Born radii");
1923     EETYPE("gb-algorithm", ir->gb_algorithm, egb_names);
1924     CTYPE ("Frequency of calculating the Born radii inside rlist");
1925     ITYPE ("nstgbradii", ir->nstgbradii, 1);
1926     CTYPE ("Cutoff for Born radii calculation; the contribution from atoms");
1927     CTYPE ("between rlist and rgbradii is updated every nstlist steps");
1928     RTYPE ("rgbradii",  ir->rgbradii, 1.0);
1929     CTYPE ("Dielectric coefficient of the implicit solvent");
1930     RTYPE ("gb-epsilon-solvent", ir->gb_epsilon_solvent, 80.0);
1931     CTYPE ("Salt concentration in M for Generalized Born models");
1932     RTYPE ("gb-saltconc",  ir->gb_saltconc, 0.0);
1933     CTYPE ("Scaling factors used in the OBC GB model. Default values are OBC(II)");
1934     RTYPE ("gb-obc-alpha", ir->gb_obc_alpha, 1.0);
1935     RTYPE ("gb-obc-beta", ir->gb_obc_beta, 0.8);
1936     RTYPE ("gb-obc-gamma", ir->gb_obc_gamma, 4.85);
1937     RTYPE ("gb-dielectric-offset", ir->gb_dielectric_offset, 0.009);
1938     EETYPE("sa-algorithm", ir->sa_algorithm, esa_names);
1939     CTYPE ("Surface tension (kJ/mol/nm^2) for the SA (nonpolar surface) part of GBSA");
1940     CTYPE ("The value -1 will set default value for Still/HCT/OBC GB-models.");
1941     RTYPE ("sa-surface-tension", ir->sa_surface_tension, -1);
1942
1943     /* Coupling stuff */
1944     CCTYPE ("OPTIONS FOR WEAK COUPLING ALGORITHMS");
1945     CTYPE ("Temperature coupling");
1946     EETYPE("tcoupl",  ir->etc,        etcoupl_names);
1947     ITYPE ("nsttcouple", ir->nsttcouple,  -1);
1948     ITYPE("nh-chain-length",     ir->opts.nhchainlength, 10);
1949     EETYPE("print-nose-hoover-chain-variables", ir->bPrintNHChains, yesno_names);
1950     CTYPE ("Groups to couple separately");
1951     STYPE ("tc-grps",     is->tcgrps,         NULL);
1952     CTYPE ("Time constant (ps) and reference temperature (K)");
1953     STYPE ("tau-t",   is->tau_t,      NULL);
1954     STYPE ("ref-t",   is->ref_t,      NULL);
1955     CTYPE ("pressure coupling");
1956     EETYPE("pcoupl",  ir->epc,        epcoupl_names);
1957     EETYPE("pcoupltype",  ir->epct,       epcoupltype_names);
1958     ITYPE ("nstpcouple", ir->nstpcouple,  -1);
1959     CTYPE ("Time constant (ps), compressibility (1/bar) and reference P (bar)");
1960     RTYPE ("tau-p",   ir->tau_p,  1.0);
1961     STYPE ("compressibility", dumstr[0],  NULL);
1962     STYPE ("ref-p",       dumstr[1],      NULL);
1963     CTYPE ("Scaling of reference coordinates, No, All or COM");
1964     EETYPE ("refcoord-scaling", ir->refcoord_scaling, erefscaling_names);
1965
1966     /* QMMM */
1967     CCTYPE ("OPTIONS FOR QMMM calculations");
1968     EETYPE("QMMM", ir->bQMMM, yesno_names);
1969     CTYPE ("Groups treated Quantum Mechanically");
1970     STYPE ("QMMM-grps",  is->QMMM,          NULL);
1971     CTYPE ("QM method");
1972     STYPE("QMmethod",     is->QMmethod, NULL);
1973     CTYPE ("QMMM scheme");
1974     EETYPE("QMMMscheme",  ir->QMMMscheme,    eQMMMscheme_names);
1975     CTYPE ("QM basisset");
1976     STYPE("QMbasis",      is->QMbasis, NULL);
1977     CTYPE ("QM charge");
1978     STYPE ("QMcharge",    is->QMcharge, NULL);
1979     CTYPE ("QM multiplicity");
1980     STYPE ("QMmult",      is->QMmult, NULL);
1981     CTYPE ("Surface Hopping");
1982     STYPE ("SH",          is->bSH, NULL);
1983     CTYPE ("CAS space options");
1984     STYPE ("CASorbitals",      is->CASorbitals,   NULL);
1985     STYPE ("CASelectrons",     is->CASelectrons,  NULL);
1986     STYPE ("SAon", is->SAon, NULL);
1987     STYPE ("SAoff", is->SAoff, NULL);
1988     STYPE ("SAsteps", is->SAsteps, NULL);
1989     CTYPE ("Scale factor for MM charges");
1990     RTYPE ("MMChargeScaleFactor", ir->scalefactor, 1.0);
1991     CTYPE ("Optimization of QM subsystem");
1992     STYPE ("bOPT",          is->bOPT, NULL);
1993     STYPE ("bTS",          is->bTS, NULL);
1994
1995     /* Simulated annealing */
1996     CCTYPE("SIMULATED ANNEALING");
1997     CTYPE ("Type of annealing for each temperature group (no/single/periodic)");
1998     STYPE ("annealing",   is->anneal,      NULL);
1999     CTYPE ("Number of time points to use for specifying annealing in each group");
2000     STYPE ("annealing-npoints", is->anneal_npoints, NULL);
2001     CTYPE ("List of times at the annealing points for each group");
2002     STYPE ("annealing-time",       is->anneal_time,       NULL);
2003     CTYPE ("Temp. at each annealing point, for each group.");
2004     STYPE ("annealing-temp",  is->anneal_temp,  NULL);
2005
2006     /* Startup run */
2007     CCTYPE ("GENERATE VELOCITIES FOR STARTUP RUN");
2008     EETYPE("gen-vel",     opts->bGenVel,  yesno_names);
2009     RTYPE ("gen-temp",    opts->tempi,    300.0);
2010     ITYPE ("gen-seed",    opts->seed,     -1);
2011
2012     /* Shake stuff */
2013     CCTYPE ("OPTIONS FOR BONDS");
2014     EETYPE("constraints", opts->nshake,   constraints);
2015     CTYPE ("Type of constraint algorithm");
2016     EETYPE("constraint-algorithm",  ir->eConstrAlg, econstr_names);
2017     CTYPE ("Do not constrain the start configuration");
2018     EETYPE("continuation", ir->bContinuation, yesno_names);
2019     CTYPE ("Use successive overrelaxation to reduce the number of shake iterations");
2020     EETYPE("Shake-SOR", ir->bShakeSOR, yesno_names);
2021     CTYPE ("Relative tolerance of shake");
2022     RTYPE ("shake-tol", ir->shake_tol, 0.0001);
2023     CTYPE ("Highest order in the expansion of the constraint coupling matrix");
2024     ITYPE ("lincs-order", ir->nProjOrder, 4);
2025     CTYPE ("Number of iterations in the final step of LINCS. 1 is fine for");
2026     CTYPE ("normal simulations, but use 2 to conserve energy in NVE runs.");
2027     CTYPE ("For energy minimization with constraints it should be 4 to 8.");
2028     ITYPE ("lincs-iter", ir->nLincsIter, 1);
2029     CTYPE ("Lincs will write a warning to the stderr if in one step a bond");
2030     CTYPE ("rotates over more degrees than");
2031     RTYPE ("lincs-warnangle", ir->LincsWarnAngle, 30.0);
2032     CTYPE ("Convert harmonic bonds to morse potentials");
2033     EETYPE("morse",       opts->bMorse, yesno_names);
2034
2035     /* Energy group exclusions */
2036     CCTYPE ("ENERGY GROUP EXCLUSIONS");
2037     CTYPE ("Pairs of energy groups for which all non-bonded interactions are excluded");
2038     STYPE ("energygrp-excl", is->egpexcl,     NULL);
2039
2040     /* Walls */
2041     CCTYPE ("WALLS");
2042     CTYPE ("Number of walls, type, atom types, densities and box-z scale factor for Ewald");
2043     ITYPE ("nwall", ir->nwall, 0);
2044     EETYPE("wall-type",     ir->wall_type,   ewt_names);
2045     RTYPE ("wall-r-linpot", ir->wall_r_linpot, -1);
2046     STYPE ("wall-atomtype", is->wall_atomtype, NULL);
2047     STYPE ("wall-density",  is->wall_density,  NULL);
2048     RTYPE ("wall-ewald-zfac", ir->wall_ewald_zfac, 3);
2049
2050     /* COM pulling */
2051     CCTYPE("COM PULLING");
2052     CTYPE("Pull type: no, umbrella, constraint or constant-force");
2053     EETYPE("pull",          ir->ePull, epull_names);
2054     if (ir->ePull != epullNO)
2055     {
2056         snew(ir->pull, 1);
2057         is->pull_grp = read_pullparams(&ninp, &inp, ir->pull, &opts->pull_start, wi);
2058     }
2059
2060     /* Enforced rotation */
2061     CCTYPE("ENFORCED ROTATION");
2062     CTYPE("Enforced rotation: No or Yes");
2063     EETYPE("rotation",       ir->bRot, yesno_names);
2064     if (ir->bRot)
2065     {
2066         snew(ir->rot, 1);
2067         is->rot_grp = read_rotparams(&ninp, &inp, ir->rot, wi);
2068     }
2069
2070     /* Interactive MD */
2071     ir->bIMD = FALSE;
2072     CCTYPE("Group to display and/or manipulate in interactive MD session");
2073     STYPE ("IMD-group", is->imd_grp, NULL);
2074     if (is->imd_grp[0] != '\0')
2075     {
2076         snew(ir->imd, 1);
2077         ir->bIMD = TRUE;
2078     }
2079
2080     /* Refinement */
2081     CCTYPE("NMR refinement stuff");
2082     CTYPE ("Distance restraints type: No, Simple or Ensemble");
2083     EETYPE("disre",       ir->eDisre,     edisre_names);
2084     CTYPE ("Force weighting of pairs in one distance restraint: Conservative or Equal");
2085     EETYPE("disre-weighting", ir->eDisreWeighting, edisreweighting_names);
2086     CTYPE ("Use sqrt of the time averaged times the instantaneous violation");
2087     EETYPE("disre-mixed", ir->bDisreMixed, yesno_names);
2088     RTYPE ("disre-fc",    ir->dr_fc,  1000.0);
2089     RTYPE ("disre-tau",   ir->dr_tau, 0.0);
2090     CTYPE ("Output frequency for pair distances to energy file");
2091     ITYPE ("nstdisreout", ir->nstdisreout, 100);
2092     CTYPE ("Orientation restraints: No or Yes");
2093     EETYPE("orire",       opts->bOrire,   yesno_names);
2094     CTYPE ("Orientation restraints force constant and tau for time averaging");
2095     RTYPE ("orire-fc",    ir->orires_fc,  0.0);
2096     RTYPE ("orire-tau",   ir->orires_tau, 0.0);
2097     STYPE ("orire-fitgrp", is->orirefitgrp,    NULL);
2098     CTYPE ("Output frequency for trace(SD) and S to energy file");
2099     ITYPE ("nstorireout", ir->nstorireout, 100);
2100
2101     /* free energy variables */
2102     CCTYPE ("Free energy variables");
2103     EETYPE("free-energy", ir->efep, efep_names);
2104     STYPE ("couple-moltype",  is->couple_moltype,  NULL);
2105     EETYPE("couple-lambda0", opts->couple_lam0, couple_lam);
2106     EETYPE("couple-lambda1", opts->couple_lam1, couple_lam);
2107     EETYPE("couple-intramol", opts->bCoupleIntra, yesno_names);
2108
2109     RTYPE ("init-lambda", fep->init_lambda, -1); /* start with -1 so
2110                                                     we can recognize if
2111                                                     it was not entered */
2112     ITYPE ("init-lambda-state", fep->init_fep_state, -1);
2113     RTYPE ("delta-lambda", fep->delta_lambda, 0.0);
2114     ITYPE ("nstdhdl", fep->nstdhdl, 50);
2115     STYPE ("fep-lambdas", is->fep_lambda[efptFEP], NULL);
2116     STYPE ("mass-lambdas", is->fep_lambda[efptMASS], NULL);
2117     STYPE ("coul-lambdas", is->fep_lambda[efptCOUL], NULL);
2118     STYPE ("vdw-lambdas", is->fep_lambda[efptVDW], NULL);
2119     STYPE ("bonded-lambdas", is->fep_lambda[efptBONDED], NULL);
2120     STYPE ("restraint-lambdas", is->fep_lambda[efptRESTRAINT], NULL);
2121     STYPE ("temperature-lambdas", is->fep_lambda[efptTEMPERATURE], NULL);
2122     ITYPE ("calc-lambda-neighbors", fep->lambda_neighbors, 1);
2123     STYPE ("init-lambda-weights", is->lambda_weights, NULL);
2124     EETYPE("dhdl-print-energy", fep->bPrintEnergy, yesno_names);
2125     RTYPE ("sc-alpha", fep->sc_alpha, 0.0);
2126     ITYPE ("sc-power", fep->sc_power, 1);
2127     RTYPE ("sc-r-power", fep->sc_r_power, 6.0);
2128     RTYPE ("sc-sigma", fep->sc_sigma, 0.3);
2129     EETYPE("sc-coul", fep->bScCoul, yesno_names);
2130     ITYPE ("dh_hist_size", fep->dh_hist_size, 0);
2131     RTYPE ("dh_hist_spacing", fep->dh_hist_spacing, 0.1);
2132     EETYPE("separate-dhdl-file", fep->separate_dhdl_file,
2133            separate_dhdl_file_names);
2134     EETYPE("dhdl-derivatives", fep->dhdl_derivatives, dhdl_derivatives_names);
2135     ITYPE ("dh_hist_size", fep->dh_hist_size, 0);
2136     RTYPE ("dh_hist_spacing", fep->dh_hist_spacing, 0.1);
2137
2138     /* Non-equilibrium MD stuff */
2139     CCTYPE("Non-equilibrium MD stuff");
2140     STYPE ("acc-grps",    is->accgrps,        NULL);
2141     STYPE ("accelerate",  is->acc,            NULL);
2142     STYPE ("freezegrps",  is->freeze,         NULL);
2143     STYPE ("freezedim",   is->frdim,          NULL);
2144     RTYPE ("cos-acceleration", ir->cos_accel, 0);
2145     STYPE ("deform",      is->deform,         NULL);
2146
2147     /* simulated tempering variables */
2148     CCTYPE("simulated tempering variables");
2149     EETYPE("simulated-tempering", ir->bSimTemp, yesno_names);
2150     EETYPE("simulated-tempering-scaling", ir->simtempvals->eSimTempScale, esimtemp_names);
2151     RTYPE("sim-temp-low", ir->simtempvals->simtemp_low, 300.0);
2152     RTYPE("sim-temp-high", ir->simtempvals->simtemp_high, 300.0);
2153
2154     /* expanded ensemble variables */
2155     if (ir->efep == efepEXPANDED || ir->bSimTemp)
2156     {
2157         read_expandedparams(&ninp, &inp, expand, wi);
2158     }
2159
2160     /* Electric fields */
2161     CCTYPE("Electric fields");
2162     CTYPE ("Format is number of terms (int) and for all terms an amplitude (real)");
2163     CTYPE ("and a phase angle (real)");
2164     STYPE ("E-x",     is->efield_x,   NULL);
2165     STYPE ("E-xt",    is->efield_xt,  NULL);
2166     STYPE ("E-y",     is->efield_y,   NULL);
2167     STYPE ("E-yt",    is->efield_yt,  NULL);
2168     STYPE ("E-z",     is->efield_z,   NULL);
2169     STYPE ("E-zt",    is->efield_zt,  NULL);
2170
2171     CCTYPE("Ion/water position swapping for computational electrophysiology setups");
2172     CTYPE("Swap positions along direction: no, X, Y, Z");
2173     EETYPE("swapcoords", ir->eSwapCoords, eSwapTypes_names);
2174     if (ir->eSwapCoords != eswapNO)
2175     {
2176         snew(ir->swap, 1);
2177         CTYPE("Swap attempt frequency");
2178         ITYPE("swap-frequency", ir->swap->nstswap, 1);
2179         CTYPE("Two index groups that contain the compartment-partitioning atoms");
2180         STYPE("split-group0", splitgrp0, NULL);
2181         STYPE("split-group1", splitgrp1, NULL);
2182         CTYPE("Use center of mass of split groups (yes/no), otherwise center of geometry is used");
2183         EETYPE("massw-split0", ir->swap->massw_split[0], yesno_names);
2184         EETYPE("massw-split1", ir->swap->massw_split[1], yesno_names);
2185
2186         CTYPE("Group name of ions that can be exchanged with solvent molecules");
2187         STYPE("swap-group", swapgrp, NULL);
2188         CTYPE("Group name of solvent molecules");
2189         STYPE("solvent-group", solgrp, NULL);
2190
2191         CTYPE("Split cylinder: radius, upper and lower extension (nm) (this will define the channels)");
2192         CTYPE("Note that the split cylinder settings do not have an influence on the swapping protocol,");
2193         CTYPE("however, if correctly defined, the ion permeation events are counted per channel");
2194         RTYPE("cyl0-r", ir->swap->cyl0r, 2.0);
2195         RTYPE("cyl0-up", ir->swap->cyl0u, 1.0);
2196         RTYPE("cyl0-down", ir->swap->cyl0l, 1.0);
2197         RTYPE("cyl1-r", ir->swap->cyl1r, 2.0);
2198         RTYPE("cyl1-up", ir->swap->cyl1u, 1.0);
2199         RTYPE("cyl1-down", ir->swap->cyl1l, 1.0);
2200
2201         CTYPE("Average the number of ions per compartment over these many swap attempt steps");
2202         ITYPE("coupl-steps", ir->swap->nAverage, 10);
2203         CTYPE("Requested number of anions and cations for each of the two compartments");
2204         CTYPE("-1 means fix the numbers as found in time step 0");
2205         ITYPE("anionsA", ir->swap->nanions[0], -1);
2206         ITYPE("cationsA", ir->swap->ncations[0], -1);
2207         ITYPE("anionsB", ir->swap->nanions[1], -1);
2208         ITYPE("cationsB", ir->swap->ncations[1], -1);
2209         CTYPE("Start to swap ions if threshold difference to requested count is reached");
2210         RTYPE("threshold", ir->swap->threshold, 1.0);
2211     }
2212
2213     /* AdResS defined thingies */
2214     CCTYPE ("AdResS parameters");
2215     EETYPE("adress",       ir->bAdress, yesno_names);
2216     if (ir->bAdress)
2217     {
2218         snew(ir->adress, 1);
2219         read_adressparams(&ninp, &inp, ir->adress, wi);
2220     }
2221
2222     /* User defined thingies */
2223     CCTYPE ("User defined thingies");
2224     STYPE ("user1-grps",  is->user1,          NULL);
2225     STYPE ("user2-grps",  is->user2,          NULL);
2226     ITYPE ("userint1",    ir->userint1,   0);
2227     ITYPE ("userint2",    ir->userint2,   0);
2228     ITYPE ("userint3",    ir->userint3,   0);
2229     ITYPE ("userint4",    ir->userint4,   0);
2230     RTYPE ("userreal1",   ir->userreal1,  0);
2231     RTYPE ("userreal2",   ir->userreal2,  0);
2232     RTYPE ("userreal3",   ir->userreal3,  0);
2233     RTYPE ("userreal4",   ir->userreal4,  0);
2234 #undef CTYPE
2235
2236     write_inpfile(mdparout, ninp, inp, FALSE, wi);
2237     for (i = 0; (i < ninp); i++)
2238     {
2239         sfree(inp[i].name);
2240         sfree(inp[i].value);
2241     }
2242     sfree(inp);
2243
2244     /* Process options if necessary */
2245     for (m = 0; m < 2; m++)
2246     {
2247         for (i = 0; i < 2*DIM; i++)
2248         {
2249             dumdub[m][i] = 0.0;
2250         }
2251         if (ir->epc)
2252         {
2253             switch (ir->epct)
2254             {
2255                 case epctISOTROPIC:
2256                     if (sscanf(dumstr[m], "%lf", &(dumdub[m][XX])) != 1)
2257                     {
2258                         warning_error(wi, "Pressure coupling not enough values (I need 1)");
2259                     }
2260                     dumdub[m][YY] = dumdub[m][ZZ] = dumdub[m][XX];
2261                     break;
2262                 case epctSEMIISOTROPIC:
2263                 case epctSURFACETENSION:
2264                     if (sscanf(dumstr[m], "%lf%lf",
2265                                &(dumdub[m][XX]), &(dumdub[m][ZZ])) != 2)
2266                     {
2267                         warning_error(wi, "Pressure coupling not enough values (I need 2)");
2268                     }
2269                     dumdub[m][YY] = dumdub[m][XX];
2270                     break;
2271                 case epctANISOTROPIC:
2272                     if (sscanf(dumstr[m], "%lf%lf%lf%lf%lf%lf",
2273                                &(dumdub[m][XX]), &(dumdub[m][YY]), &(dumdub[m][ZZ]),
2274                                &(dumdub[m][3]), &(dumdub[m][4]), &(dumdub[m][5])) != 6)
2275                     {
2276                         warning_error(wi, "Pressure coupling not enough values (I need 6)");
2277                     }
2278                     break;
2279                 default:
2280                     gmx_fatal(FARGS, "Pressure coupling type %s not implemented yet",
2281                               epcoupltype_names[ir->epct]);
2282             }
2283         }
2284     }
2285     clear_mat(ir->ref_p);
2286     clear_mat(ir->compress);
2287     for (i = 0; i < DIM; i++)
2288     {
2289         ir->ref_p[i][i]    = dumdub[1][i];
2290         ir->compress[i][i] = dumdub[0][i];
2291     }
2292     if (ir->epct == epctANISOTROPIC)
2293     {
2294         ir->ref_p[XX][YY] = dumdub[1][3];
2295         ir->ref_p[XX][ZZ] = dumdub[1][4];
2296         ir->ref_p[YY][ZZ] = dumdub[1][5];
2297         if (ir->ref_p[XX][YY] != 0 && ir->ref_p[XX][ZZ] != 0 && ir->ref_p[YY][ZZ] != 0)
2298         {
2299             warning(wi, "All off-diagonal reference pressures are non-zero. Are you sure you want to apply a threefold shear stress?\n");
2300         }
2301         ir->compress[XX][YY] = dumdub[0][3];
2302         ir->compress[XX][ZZ] = dumdub[0][4];
2303         ir->compress[YY][ZZ] = dumdub[0][5];
2304         for (i = 0; i < DIM; i++)
2305         {
2306             for (m = 0; m < i; m++)
2307             {
2308                 ir->ref_p[i][m]    = ir->ref_p[m][i];
2309                 ir->compress[i][m] = ir->compress[m][i];
2310             }
2311         }
2312     }
2313
2314     if (ir->comm_mode == ecmNO)
2315     {
2316         ir->nstcomm = 0;
2317     }
2318
2319     opts->couple_moltype = NULL;
2320     if (strlen(is->couple_moltype) > 0)
2321     {
2322         if (ir->efep != efepNO)
2323         {
2324             opts->couple_moltype = strdup(is->couple_moltype);
2325             if (opts->couple_lam0 == opts->couple_lam1)
2326             {
2327                 warning(wi, "The lambda=0 and lambda=1 states for coupling are identical");
2328             }
2329             if (ir->eI == eiMD && (opts->couple_lam0 == ecouplamNONE ||
2330                                    opts->couple_lam1 == ecouplamNONE))
2331             {
2332                 warning(wi, "For proper sampling of the (nearly) decoupled state, stochastic dynamics should be used");
2333             }
2334         }
2335         else
2336         {
2337             warning(wi, "Can not couple a molecule with free_energy = no");
2338         }
2339     }
2340     /* FREE ENERGY AND EXPANDED ENSEMBLE OPTIONS */
2341     if (ir->efep != efepNO)
2342     {
2343         if (fep->delta_lambda > 0)
2344         {
2345             ir->efep = efepSLOWGROWTH;
2346         }
2347     }
2348
2349     if (ir->bSimTemp)
2350     {
2351         fep->bPrintEnergy = TRUE;
2352         /* always print out the energy to dhdl if we are doing expanded ensemble, since we need the total energy
2353            if the temperature is changing. */
2354     }
2355
2356     if ((ir->efep != efepNO) || ir->bSimTemp)
2357     {
2358         ir->bExpanded = FALSE;
2359         if ((ir->efep == efepEXPANDED) || ir->bSimTemp)
2360         {
2361             ir->bExpanded = TRUE;
2362         }
2363         do_fep_params(ir, is->fep_lambda, is->lambda_weights);
2364         if (ir->bSimTemp) /* done after fep params */
2365         {
2366             do_simtemp_params(ir);
2367         }
2368     }
2369     else
2370     {
2371         ir->fepvals->n_lambda = 0;
2372     }
2373
2374     /* WALL PARAMETERS */
2375
2376     do_wall_params(ir, is->wall_atomtype, is->wall_density, opts);
2377
2378     /* ORIENTATION RESTRAINT PARAMETERS */
2379
2380     if (opts->bOrire && str_nelem(is->orirefitgrp, MAXPTR, NULL) != 1)
2381     {
2382         warning_error(wi, "ERROR: Need one orientation restraint fit group\n");
2383     }
2384
2385     /* DEFORMATION PARAMETERS */
2386
2387     clear_mat(ir->deform);
2388     for (i = 0; i < 6; i++)
2389     {
2390         dumdub[0][i] = 0;
2391     }
2392     m = sscanf(is->deform, "%lf %lf %lf %lf %lf %lf",
2393                &(dumdub[0][0]), &(dumdub[0][1]), &(dumdub[0][2]),
2394                &(dumdub[0][3]), &(dumdub[0][4]), &(dumdub[0][5]));
2395     for (i = 0; i < 3; i++)
2396     {
2397         ir->deform[i][i] = dumdub[0][i];
2398     }
2399     ir->deform[YY][XX] = dumdub[0][3];
2400     ir->deform[ZZ][XX] = dumdub[0][4];
2401     ir->deform[ZZ][YY] = dumdub[0][5];
2402     if (ir->epc != epcNO)
2403     {
2404         for (i = 0; i < 3; i++)
2405         {
2406             for (j = 0; j <= i; j++)
2407             {
2408                 if (ir->deform[i][j] != 0 && ir->compress[i][j] != 0)
2409                 {
2410                     warning_error(wi, "A box element has deform set and compressibility > 0");
2411                 }
2412             }
2413         }
2414         for (i = 0; i < 3; i++)
2415         {
2416             for (j = 0; j < i; j++)
2417             {
2418                 if (ir->deform[i][j] != 0)
2419                 {
2420                     for (m = j; m < DIM; m++)
2421                     {
2422                         if (ir->compress[m][j] != 0)
2423                         {
2424                             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.");
2425                             warning(wi, warn_buf);
2426                         }
2427                     }
2428                 }
2429             }
2430         }
2431     }
2432
2433     /* Ion/water position swapping checks */
2434     if (ir->eSwapCoords != eswapNO)
2435     {
2436         if (ir->swap->nstswap < 1)
2437         {
2438             warning_error(wi, "swap_frequency must be 1 or larger when ion swapping is requested");
2439         }
2440         if (ir->swap->nAverage < 1)
2441         {
2442             warning_error(wi, "coupl_steps must be 1 or larger.\n");
2443         }
2444         if (ir->swap->threshold < 1.0)
2445         {
2446             warning_error(wi, "Ion count threshold must be at least 1.\n");
2447         }
2448     }
2449
2450     sfree(dumstr[0]);
2451     sfree(dumstr[1]);
2452 }
2453
2454 static int search_QMstring(const char *s, int ng, const char *gn[])
2455 {
2456     /* same as normal search_string, but this one searches QM strings */
2457     int i;
2458
2459     for (i = 0; (i < ng); i++)
2460     {
2461         if (gmx_strcasecmp(s, gn[i]) == 0)
2462         {
2463             return i;
2464         }
2465     }
2466
2467     gmx_fatal(FARGS, "this QM method or basisset (%s) is not implemented\n!", s);
2468
2469     return -1;
2470
2471 } /* search_QMstring */
2472
2473 /* We would like gn to be const as well, but C doesn't allow this */
2474 int search_string(const char *s, int ng, char *gn[])
2475 {
2476     int i;
2477
2478     for (i = 0; (i < ng); i++)
2479     {
2480         if (gmx_strcasecmp(s, gn[i]) == 0)
2481         {
2482             return i;
2483         }
2484     }
2485
2486     gmx_fatal(FARGS,
2487               "Group %s referenced in the .mdp file was not found in the index file.\n"
2488               "Group names must match either [moleculetype] names or custom index group\n"
2489               "names, in which case you must supply an index file to the '-n' option\n"
2490               "of grompp.",
2491               s);
2492
2493     return -1;
2494 }
2495
2496 static gmx_bool do_numbering(int natoms, gmx_groups_t *groups, int ng, char *ptrs[],
2497                              t_blocka *block, char *gnames[],
2498                              int gtype, int restnm,
2499                              int grptp, gmx_bool bVerbose,
2500                              warninp_t wi)
2501 {
2502     unsigned short *cbuf;
2503     t_grps         *grps = &(groups->grps[gtype]);
2504     int             i, j, gid, aj, ognr, ntot = 0;
2505     const char     *title;
2506     gmx_bool        bRest;
2507     char            warn_buf[STRLEN];
2508
2509     if (debug)
2510     {
2511         fprintf(debug, "Starting numbering %d groups of type %d\n", ng, gtype);
2512     }
2513
2514     title = gtypes[gtype];
2515
2516     snew(cbuf, natoms);
2517     /* Mark all id's as not set */
2518     for (i = 0; (i < natoms); i++)
2519     {
2520         cbuf[i] = NOGID;
2521     }
2522
2523     snew(grps->nm_ind, ng+1); /* +1 for possible rest group */
2524     for (i = 0; (i < ng); i++)
2525     {
2526         /* Lookup the group name in the block structure */
2527         gid = search_string(ptrs[i], block->nr, gnames);
2528         if ((grptp != egrptpONE) || (i == 0))
2529         {
2530             grps->nm_ind[grps->nr++] = gid;
2531         }
2532         if (debug)
2533         {
2534             fprintf(debug, "Found gid %d for group %s\n", gid, ptrs[i]);
2535         }
2536
2537         /* Now go over the atoms in the group */
2538         for (j = block->index[gid]; (j < block->index[gid+1]); j++)
2539         {
2540
2541             aj = block->a[j];
2542
2543             /* Range checking */
2544             if ((aj < 0) || (aj >= natoms))
2545             {
2546                 gmx_fatal(FARGS, "Invalid atom number %d in indexfile", aj);
2547             }
2548             /* Lookup up the old group number */
2549             ognr = cbuf[aj];
2550             if (ognr != NOGID)
2551             {
2552                 gmx_fatal(FARGS, "Atom %d in multiple %s groups (%d and %d)",
2553                           aj+1, title, ognr+1, i+1);
2554             }
2555             else
2556             {
2557                 /* Store the group number in buffer */
2558                 if (grptp == egrptpONE)
2559                 {
2560                     cbuf[aj] = 0;
2561                 }
2562                 else
2563                 {
2564                     cbuf[aj] = i;
2565                 }
2566                 ntot++;
2567             }
2568         }
2569     }
2570
2571     /* Now check whether we have done all atoms */
2572     bRest = FALSE;
2573     if (ntot != natoms)
2574     {
2575         if (grptp == egrptpALL)
2576         {
2577             gmx_fatal(FARGS, "%d atoms are not part of any of the %s groups",
2578                       natoms-ntot, title);
2579         }
2580         else if (grptp == egrptpPART)
2581         {
2582             sprintf(warn_buf, "%d atoms are not part of any of the %s groups",
2583                     natoms-ntot, title);
2584             warning_note(wi, warn_buf);
2585         }
2586         /* Assign all atoms currently unassigned to a rest group */
2587         for (j = 0; (j < natoms); j++)
2588         {
2589             if (cbuf[j] == NOGID)
2590             {
2591                 cbuf[j] = grps->nr;
2592                 bRest   = TRUE;
2593             }
2594         }
2595         if (grptp != egrptpPART)
2596         {
2597             if (bVerbose)
2598             {
2599                 fprintf(stderr,
2600                         "Making dummy/rest group for %s containing %d elements\n",
2601                         title, natoms-ntot);
2602             }
2603             /* Add group name "rest" */
2604             grps->nm_ind[grps->nr] = restnm;
2605
2606             /* Assign the rest name to all atoms not currently assigned to a group */
2607             for (j = 0; (j < natoms); j++)
2608             {
2609                 if (cbuf[j] == NOGID)
2610                 {
2611                     cbuf[j] = grps->nr;
2612                 }
2613             }
2614             grps->nr++;
2615         }
2616     }
2617
2618     if (grps->nr == 1 && (ntot == 0 || ntot == natoms))
2619     {
2620         /* All atoms are part of one (or no) group, no index required */
2621         groups->ngrpnr[gtype] = 0;
2622         groups->grpnr[gtype]  = NULL;
2623     }
2624     else
2625     {
2626         groups->ngrpnr[gtype] = natoms;
2627         snew(groups->grpnr[gtype], natoms);
2628         for (j = 0; (j < natoms); j++)
2629         {
2630             groups->grpnr[gtype][j] = cbuf[j];
2631         }
2632     }
2633
2634     sfree(cbuf);
2635
2636     return (bRest && grptp == egrptpPART);
2637 }
2638
2639 static void calc_nrdf(gmx_mtop_t *mtop, t_inputrec *ir, char **gnames)
2640 {
2641     t_grpopts              *opts;
2642     gmx_groups_t           *groups;
2643     t_pull                 *pull;
2644     int                     natoms, ai, aj, i, j, d, g, imin, jmin;
2645     t_iatom                *ia;
2646     int                    *nrdf2, *na_vcm, na_tot;
2647     double                 *nrdf_tc, *nrdf_vcm, nrdf_uc, n_sub = 0;
2648     gmx_mtop_atomloop_all_t aloop;
2649     t_atom                 *atom;
2650     int                     mb, mol, ftype, as;
2651     gmx_molblock_t         *molb;
2652     gmx_moltype_t          *molt;
2653
2654     /* Calculate nrdf.
2655      * First calc 3xnr-atoms for each group
2656      * then subtract half a degree of freedom for each constraint
2657      *
2658      * Only atoms and nuclei contribute to the degrees of freedom...
2659      */
2660
2661     opts = &ir->opts;
2662
2663     groups = &mtop->groups;
2664     natoms = mtop->natoms;
2665
2666     /* Allocate one more for a possible rest group */
2667     /* We need to sum degrees of freedom into doubles,
2668      * since floats give too low nrdf's above 3 million atoms.
2669      */
2670     snew(nrdf_tc, groups->grps[egcTC].nr+1);
2671     snew(nrdf_vcm, groups->grps[egcVCM].nr+1);
2672     snew(na_vcm, groups->grps[egcVCM].nr+1);
2673
2674     for (i = 0; i < groups->grps[egcTC].nr; i++)
2675     {
2676         nrdf_tc[i] = 0;
2677     }
2678     for (i = 0; i < groups->grps[egcVCM].nr+1; i++)
2679     {
2680         nrdf_vcm[i] = 0;
2681     }
2682
2683     snew(nrdf2, natoms);
2684     aloop = gmx_mtop_atomloop_all_init(mtop);
2685     while (gmx_mtop_atomloop_all_next(aloop, &i, &atom))
2686     {
2687         nrdf2[i] = 0;
2688         if (atom->ptype == eptAtom || atom->ptype == eptNucleus)
2689         {
2690             g = ggrpnr(groups, egcFREEZE, i);
2691             /* Double count nrdf for particle i */
2692             for (d = 0; d < DIM; d++)
2693             {
2694                 if (opts->nFreeze[g][d] == 0)
2695                 {
2696                     nrdf2[i] += 2;
2697                 }
2698             }
2699             nrdf_tc [ggrpnr(groups, egcTC, i)]  += 0.5*nrdf2[i];
2700             nrdf_vcm[ggrpnr(groups, egcVCM, i)] += 0.5*nrdf2[i];
2701         }
2702     }
2703
2704     as = 0;
2705     for (mb = 0; mb < mtop->nmolblock; mb++)
2706     {
2707         molb = &mtop->molblock[mb];
2708         molt = &mtop->moltype[molb->type];
2709         atom = molt->atoms.atom;
2710         for (mol = 0; mol < molb->nmol; mol++)
2711         {
2712             for (ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++)
2713             {
2714                 ia = molt->ilist[ftype].iatoms;
2715                 for (i = 0; i < molt->ilist[ftype].nr; )
2716                 {
2717                     /* Subtract degrees of freedom for the constraints,
2718                      * if the particles still have degrees of freedom left.
2719                      * If one of the particles is a vsite or a shell, then all
2720                      * constraint motion will go there, but since they do not
2721                      * contribute to the constraints the degrees of freedom do not
2722                      * change.
2723                      */
2724                     ai = as + ia[1];
2725                     aj = as + ia[2];
2726                     if (((atom[ia[1]].ptype == eptNucleus) ||
2727                          (atom[ia[1]].ptype == eptAtom)) &&
2728                         ((atom[ia[2]].ptype == eptNucleus) ||
2729                          (atom[ia[2]].ptype == eptAtom)))
2730                     {
2731                         if (nrdf2[ai] > 0)
2732                         {
2733                             jmin = 1;
2734                         }
2735                         else
2736                         {
2737                             jmin = 2;
2738                         }
2739                         if (nrdf2[aj] > 0)
2740                         {
2741                             imin = 1;
2742                         }
2743                         else
2744                         {
2745                             imin = 2;
2746                         }
2747                         imin       = min(imin, nrdf2[ai]);
2748                         jmin       = min(jmin, nrdf2[aj]);
2749                         nrdf2[ai] -= imin;
2750                         nrdf2[aj] -= jmin;
2751                         nrdf_tc [ggrpnr(groups, egcTC, ai)]  -= 0.5*imin;
2752                         nrdf_tc [ggrpnr(groups, egcTC, aj)]  -= 0.5*jmin;
2753                         nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5*imin;
2754                         nrdf_vcm[ggrpnr(groups, egcVCM, aj)] -= 0.5*jmin;
2755                     }
2756                     ia += interaction_function[ftype].nratoms+1;
2757                     i  += interaction_function[ftype].nratoms+1;
2758                 }
2759             }
2760             ia = molt->ilist[F_SETTLE].iatoms;
2761             for (i = 0; i < molt->ilist[F_SETTLE].nr; )
2762             {
2763                 /* Subtract 1 dof from every atom in the SETTLE */
2764                 for (j = 0; j < 3; j++)
2765                 {
2766                     ai         = as + ia[1+j];
2767                     imin       = min(2, nrdf2[ai]);
2768                     nrdf2[ai] -= imin;
2769                     nrdf_tc [ggrpnr(groups, egcTC, ai)]  -= 0.5*imin;
2770                     nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5*imin;
2771                 }
2772                 ia += 4;
2773                 i  += 4;
2774             }
2775             as += molt->atoms.nr;
2776         }
2777     }
2778
2779     if (ir->ePull == epullCONSTRAINT)
2780     {
2781         /* Correct nrdf for the COM constraints.
2782          * We correct using the TC and VCM group of the first atom
2783          * in the reference and pull group. If atoms in one pull group
2784          * belong to different TC or VCM groups it is anyhow difficult
2785          * to determine the optimal nrdf assignment.
2786          */
2787         pull = ir->pull;
2788
2789         for (i = 0; i < pull->ncoord; i++)
2790         {
2791             imin = 1;
2792
2793             for (j = 0; j < 2; j++)
2794             {
2795                 const t_pull_group *pgrp;
2796
2797                 pgrp = &pull->group[pull->coord[i].group[j]];
2798
2799                 if (pgrp->nat > 0)
2800                 {
2801                     /* Subtract 1/2 dof from each group */
2802                     ai = pgrp->ind[0];
2803                     nrdf_tc [ggrpnr(groups, egcTC, ai)]  -= 0.5*imin;
2804                     nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5*imin;
2805                     if (nrdf_tc[ggrpnr(groups, egcTC, ai)] < 0)
2806                     {
2807                         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)]]);
2808                     }
2809                 }
2810                 else
2811                 {
2812                     /* We need to subtract the whole DOF from group j=1 */
2813                     imin += 1;
2814                 }
2815             }
2816         }
2817     }
2818
2819     if (ir->nstcomm != 0)
2820     {
2821         /* Subtract 3 from the number of degrees of freedom in each vcm group
2822          * when com translation is removed and 6 when rotation is removed
2823          * as well.
2824          */
2825         switch (ir->comm_mode)
2826         {
2827             case ecmLINEAR:
2828                 n_sub = ndof_com(ir);
2829                 break;
2830             case ecmANGULAR:
2831                 n_sub = 6;
2832                 break;
2833             default:
2834                 n_sub = 0;
2835                 gmx_incons("Checking comm_mode");
2836         }
2837
2838         for (i = 0; i < groups->grps[egcTC].nr; i++)
2839         {
2840             /* Count the number of atoms of TC group i for every VCM group */
2841             for (j = 0; j < groups->grps[egcVCM].nr+1; j++)
2842             {
2843                 na_vcm[j] = 0;
2844             }
2845             na_tot = 0;
2846             for (ai = 0; ai < natoms; ai++)
2847             {
2848                 if (ggrpnr(groups, egcTC, ai) == i)
2849                 {
2850                     na_vcm[ggrpnr(groups, egcVCM, ai)]++;
2851                     na_tot++;
2852                 }
2853             }
2854             /* Correct for VCM removal according to the fraction of each VCM
2855              * group present in this TC group.
2856              */
2857             nrdf_uc = nrdf_tc[i];
2858             if (debug)
2859             {
2860                 fprintf(debug, "T-group[%d] nrdf_uc = %g, n_sub = %g\n",
2861                         i, nrdf_uc, n_sub);
2862             }
2863             nrdf_tc[i] = 0;
2864             for (j = 0; j < groups->grps[egcVCM].nr+1; j++)
2865             {
2866                 if (nrdf_vcm[j] > n_sub)
2867                 {
2868                     nrdf_tc[i] += nrdf_uc*((double)na_vcm[j]/(double)na_tot)*
2869                         (nrdf_vcm[j] - n_sub)/nrdf_vcm[j];
2870                 }
2871                 if (debug)
2872                 {
2873                     fprintf(debug, "  nrdf_vcm[%d] = %g, nrdf = %g\n",
2874                             j, nrdf_vcm[j], nrdf_tc[i]);
2875                 }
2876             }
2877         }
2878     }
2879     for (i = 0; (i < groups->grps[egcTC].nr); i++)
2880     {
2881         opts->nrdf[i] = nrdf_tc[i];
2882         if (opts->nrdf[i] < 0)
2883         {
2884             opts->nrdf[i] = 0;
2885         }
2886         fprintf(stderr,
2887                 "Number of degrees of freedom in T-Coupling group %s is %.2f\n",
2888                 gnames[groups->grps[egcTC].nm_ind[i]], opts->nrdf[i]);
2889     }
2890
2891     sfree(nrdf2);
2892     sfree(nrdf_tc);
2893     sfree(nrdf_vcm);
2894     sfree(na_vcm);
2895 }
2896
2897 static void decode_cos(char *s, t_cosines *cosine)
2898 {
2899     char   *t;
2900     char    format[STRLEN], f1[STRLEN];
2901     double  a, phi;
2902     int     i;
2903
2904     t = strdup(s);
2905     trim(t);
2906
2907     cosine->n   = 0;
2908     cosine->a   = NULL;
2909     cosine->phi = NULL;
2910     if (strlen(t))
2911     {
2912         sscanf(t, "%d", &(cosine->n));
2913         if (cosine->n <= 0)
2914         {
2915             cosine->n = 0;
2916         }
2917         else
2918         {
2919             snew(cosine->a, cosine->n);
2920             snew(cosine->phi, cosine->n);
2921
2922             sprintf(format, "%%*d");
2923             for (i = 0; (i < cosine->n); i++)
2924             {
2925                 strcpy(f1, format);
2926                 strcat(f1, "%lf%lf");
2927                 if (sscanf(t, f1, &a, &phi) < 2)
2928                 {
2929                     gmx_fatal(FARGS, "Invalid input for electric field shift: '%s'", t);
2930                 }
2931                 cosine->a[i]   = a;
2932                 cosine->phi[i] = phi;
2933                 strcat(format, "%*lf%*lf");
2934             }
2935         }
2936     }
2937     sfree(t);
2938 }
2939
2940 static gmx_bool do_egp_flag(t_inputrec *ir, gmx_groups_t *groups,
2941                             const char *option, const char *val, int flag)
2942 {
2943     /* The maximum number of energy group pairs would be MAXPTR*(MAXPTR+1)/2.
2944      * But since this is much larger than STRLEN, such a line can not be parsed.
2945      * The real maximum is the number of names that fit in a string: STRLEN/2.
2946      */
2947 #define EGP_MAX (STRLEN/2)
2948     int      nelem, i, j, k, nr;
2949     char    *names[EGP_MAX];
2950     char  ***gnames;
2951     gmx_bool bSet;
2952
2953     gnames = groups->grpname;
2954
2955     nelem = str_nelem(val, EGP_MAX, names);
2956     if (nelem % 2 != 0)
2957     {
2958         gmx_fatal(FARGS, "The number of groups for %s is odd", option);
2959     }
2960     nr   = groups->grps[egcENER].nr;
2961     bSet = FALSE;
2962     for (i = 0; i < nelem/2; i++)
2963     {
2964         j = 0;
2965         while ((j < nr) &&
2966                gmx_strcasecmp(names[2*i], *(gnames[groups->grps[egcENER].nm_ind[j]])))
2967         {
2968             j++;
2969         }
2970         if (j == nr)
2971         {
2972             gmx_fatal(FARGS, "%s in %s is not an energy group\n",
2973                       names[2*i], option);
2974         }
2975         k = 0;
2976         while ((k < nr) &&
2977                gmx_strcasecmp(names[2*i+1], *(gnames[groups->grps[egcENER].nm_ind[k]])))
2978         {
2979             k++;
2980         }
2981         if (k == nr)
2982         {
2983             gmx_fatal(FARGS, "%s in %s is not an energy group\n",
2984                       names[2*i+1], option);
2985         }
2986         if ((j < nr) && (k < nr))
2987         {
2988             ir->opts.egp_flags[nr*j+k] |= flag;
2989             ir->opts.egp_flags[nr*k+j] |= flag;
2990             bSet = TRUE;
2991         }
2992     }
2993
2994     return bSet;
2995 }
2996
2997
2998 static void make_swap_groups(
2999         t_swapcoords *swap,
3000         char         *swapgname,
3001         char         *splitg0name,
3002         char         *splitg1name,
3003         char         *solgname,
3004         t_blocka     *grps,
3005         char        **gnames)
3006 {
3007     int   ig = -1, i = 0, j;
3008     char *splitg;
3009
3010
3011     /* Just a quick check here, more thorough checks are in mdrun */
3012     if (strcmp(splitg0name, splitg1name) == 0)
3013     {
3014         gmx_fatal(FARGS, "The split groups can not both be '%s'.", splitg0name);
3015     }
3016
3017     /* First get the swap group index atoms */
3018     ig        = search_string(swapgname, grps->nr, gnames);
3019     swap->nat = grps->index[ig+1] - grps->index[ig];
3020     if (swap->nat > 0)
3021     {
3022         fprintf(stderr, "Swap group '%s' contains %d atoms.\n", swapgname, swap->nat);
3023         snew(swap->ind, swap->nat);
3024         for (i = 0; i < swap->nat; i++)
3025         {
3026             swap->ind[i] = grps->a[grps->index[ig]+i];
3027         }
3028     }
3029     else
3030     {
3031         gmx_fatal(FARGS, "You defined an empty group of atoms for swapping.");
3032     }
3033
3034     /* Now do so for the split groups */
3035     for (j = 0; j < 2; j++)
3036     {
3037         if (j == 0)
3038         {
3039             splitg = splitg0name;
3040         }
3041         else
3042         {
3043             splitg = splitg1name;
3044         }
3045
3046         ig                 = search_string(splitg, grps->nr, gnames);
3047         swap->nat_split[j] = grps->index[ig+1] - grps->index[ig];
3048         if (swap->nat_split[j] > 0)
3049         {
3050             fprintf(stderr, "Split group %d '%s' contains %d atom%s.\n",
3051                     j, splitg, swap->nat_split[j], (swap->nat_split[j] > 1) ? "s" : "");
3052             snew(swap->ind_split[j], swap->nat_split[j]);
3053             for (i = 0; i < swap->nat_split[j]; i++)
3054             {
3055                 swap->ind_split[j][i] = grps->a[grps->index[ig]+i];
3056             }
3057         }
3058         else
3059         {
3060             gmx_fatal(FARGS, "Split group %d has to contain at least 1 atom!", j);
3061         }
3062     }
3063
3064     /* Now get the solvent group index atoms */
3065     ig            = search_string(solgname, grps->nr, gnames);
3066     swap->nat_sol = grps->index[ig+1] - grps->index[ig];
3067     if (swap->nat_sol > 0)
3068     {
3069         fprintf(stderr, "Solvent group '%s' contains %d atoms.\n", solgname, swap->nat_sol);
3070         snew(swap->ind_sol, swap->nat_sol);
3071         for (i = 0; i < swap->nat_sol; i++)
3072         {
3073             swap->ind_sol[i] = grps->a[grps->index[ig]+i];
3074         }
3075     }
3076     else
3077     {
3078         gmx_fatal(FARGS, "You defined an empty group of solvent. Cannot exchange ions.");
3079     }
3080 }
3081
3082
3083 void make_IMD_group(t_IMD *IMDgroup, char *IMDgname, t_blocka *grps, char **gnames)
3084 {
3085     int      ig = -1, i;
3086
3087
3088     ig            = search_string(IMDgname, grps->nr, gnames);
3089     IMDgroup->nat = grps->index[ig+1] - grps->index[ig];
3090
3091     if (IMDgroup->nat > 0)
3092     {
3093         fprintf(stderr, "Group '%s' with %d atoms can be activated for interactive molecular dynamics (IMD).\n",
3094                 IMDgname, IMDgroup->nat);
3095         snew(IMDgroup->ind, IMDgroup->nat);
3096         for (i = 0; i < IMDgroup->nat; i++)
3097         {
3098             IMDgroup->ind[i] = grps->a[grps->index[ig]+i];
3099         }
3100     }
3101 }
3102
3103
3104 void do_index(const char* mdparin, const char *ndx,
3105               gmx_mtop_t *mtop,
3106               gmx_bool bVerbose,
3107               t_inputrec *ir, rvec *v,
3108               warninp_t wi)
3109 {
3110     t_blocka     *grps;
3111     gmx_groups_t *groups;
3112     int           natoms;
3113     t_symtab     *symtab;
3114     t_atoms       atoms_all;
3115     char          warnbuf[STRLEN], **gnames;
3116     int           nr, ntcg, ntau_t, nref_t, nacc, nofg, nSA, nSA_points, nSA_time, nSA_temp;
3117     real          tau_min;
3118     int           nstcmin;
3119     int           nacg, nfreeze, nfrdim, nenergy, nvcm, nuser;
3120     char         *ptr1[MAXPTR], *ptr2[MAXPTR], *ptr3[MAXPTR];
3121     int           i, j, k, restnm;
3122     real          SAtime;
3123     gmx_bool      bExcl, bTable, bSetTCpar, bAnneal, bRest;
3124     int           nQMmethod, nQMbasis, nQMcharge, nQMmult, nbSH, nCASorb, nCASelec,
3125                   nSAon, nSAoff, nSAsteps, nQMg, nbOPT, nbTS;
3126     char          warn_buf[STRLEN];
3127
3128     if (bVerbose)
3129     {
3130         fprintf(stderr, "processing index file...\n");
3131     }
3132     debug_gmx();
3133     if (ndx == NULL)
3134     {
3135         snew(grps, 1);
3136         snew(grps->index, 1);
3137         snew(gnames, 1);
3138         atoms_all = gmx_mtop_global_atoms(mtop);
3139         analyse(&atoms_all, grps, &gnames, FALSE, TRUE);
3140         free_t_atoms(&atoms_all, FALSE);
3141     }
3142     else
3143     {
3144         grps = init_index(ndx, &gnames);
3145     }
3146
3147     groups = &mtop->groups;
3148     natoms = mtop->natoms;
3149     symtab = &mtop->symtab;
3150
3151     snew(groups->grpname, grps->nr+1);
3152
3153     for (i = 0; (i < grps->nr); i++)
3154     {
3155         groups->grpname[i] = put_symtab(symtab, gnames[i]);
3156     }
3157     groups->grpname[i] = put_symtab(symtab, "rest");
3158     restnm             = i;
3159     srenew(gnames, grps->nr+1);
3160     gnames[restnm]   = *(groups->grpname[i]);
3161     groups->ngrpname = grps->nr+1;
3162
3163     set_warning_line(wi, mdparin, -1);
3164
3165     ntau_t = str_nelem(is->tau_t, MAXPTR, ptr1);
3166     nref_t = str_nelem(is->ref_t, MAXPTR, ptr2);
3167     ntcg   = str_nelem(is->tcgrps, MAXPTR, ptr3);
3168     if ((ntau_t != ntcg) || (nref_t != ntcg))
3169     {
3170         gmx_fatal(FARGS, "Invalid T coupling input: %d groups, %d ref-t values and "
3171                   "%d tau-t values", ntcg, nref_t, ntau_t);
3172     }
3173
3174     bSetTCpar = (ir->etc || EI_SD(ir->eI) || ir->eI == eiBD || EI_TPI(ir->eI));
3175     do_numbering(natoms, groups, ntcg, ptr3, grps, gnames, egcTC,
3176                  restnm, bSetTCpar ? egrptpALL : egrptpALL_GENREST, bVerbose, wi);
3177     nr            = groups->grps[egcTC].nr;
3178     ir->opts.ngtc = nr;
3179     snew(ir->opts.nrdf, nr);
3180     snew(ir->opts.tau_t, nr);
3181     snew(ir->opts.ref_t, nr);
3182     if (ir->eI == eiBD && ir->bd_fric == 0)
3183     {
3184         fprintf(stderr, "bd-fric=0, so tau-t will be used as the inverse friction constant(s)\n");
3185     }
3186
3187     if (bSetTCpar)
3188     {
3189         if (nr != nref_t)
3190         {
3191             gmx_fatal(FARGS, "Not enough ref-t and tau-t values!");
3192         }
3193
3194         tau_min = 1e20;
3195         for (i = 0; (i < nr); i++)
3196         {
3197             ir->opts.tau_t[i] = strtod(ptr1[i], NULL);
3198             if ((ir->eI == eiBD || ir->eI == eiSD2) && ir->opts.tau_t[i] <= 0)
3199             {
3200                 sprintf(warn_buf, "With integrator %s tau-t should be larger than 0", ei_names[ir->eI]);
3201                 warning_error(wi, warn_buf);
3202             }
3203
3204             if (ir->etc != etcVRESCALE && ir->opts.tau_t[i] == 0)
3205             {
3206                 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.");
3207             }
3208
3209             if (ir->opts.tau_t[i] >= 0)
3210             {
3211                 tau_min = min(tau_min, ir->opts.tau_t[i]);
3212             }
3213         }
3214         if (ir->etc != etcNO && ir->nsttcouple == -1)
3215         {
3216             ir->nsttcouple = ir_optimal_nsttcouple(ir);
3217         }
3218
3219         if (EI_VV(ir->eI))
3220         {
3221             if ((ir->etc == etcNOSEHOOVER) && (ir->epc == epcBERENDSEN))
3222             {
3223                 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");
3224             }
3225             if ((ir->epc == epcMTTK) && (ir->etc > etcNO))
3226             {
3227                 if (ir->nstpcouple != ir->nsttcouple)
3228                 {
3229                     int mincouple = min(ir->nstpcouple, ir->nsttcouple);
3230                     ir->nstpcouple = ir->nsttcouple = mincouple;
3231                     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);
3232                     warning_note(wi, warn_buf);
3233                 }
3234             }
3235         }
3236         /* velocity verlet with averaged kinetic energy KE = 0.5*(v(t+1/2) - v(t-1/2)) is implemented
3237            primarily for testing purposes, and does not work with temperature coupling other than 1 */
3238
3239         if (ETC_ANDERSEN(ir->etc))
3240         {
3241             if (ir->nsttcouple != 1)
3242             {
3243                 ir->nsttcouple = 1;
3244                 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");
3245                 warning_note(wi, warn_buf);
3246             }
3247         }
3248         nstcmin = tcouple_min_integration_steps(ir->etc);
3249         if (nstcmin > 1)
3250         {
3251             if (tau_min/(ir->delta_t*ir->nsttcouple) < nstcmin)
3252             {
3253                 sprintf(warn_buf, "For proper integration of the %s thermostat, tau-t (%g) should be at least %d times larger than nsttcouple*dt (%g)",
3254                         ETCOUPLTYPE(ir->etc),
3255                         tau_min, nstcmin,
3256                         ir->nsttcouple*ir->delta_t);
3257                 warning(wi, warn_buf);
3258             }
3259         }
3260         for (i = 0; (i < nr); i++)
3261         {
3262             ir->opts.ref_t[i] = strtod(ptr2[i], NULL);
3263             if (ir->opts.ref_t[i] < 0)
3264             {
3265                 gmx_fatal(FARGS, "ref-t for group %d negative", i);
3266             }
3267         }
3268         /* set the lambda mc temperature to the md integrator temperature (which should be defined
3269            if we are in this conditional) if mc_temp is negative */
3270         if (ir->expandedvals->mc_temp < 0)
3271         {
3272             ir->expandedvals->mc_temp = ir->opts.ref_t[0]; /*for now, set to the first reft */
3273         }
3274     }
3275
3276     /* Simulated annealing for each group. There are nr groups */
3277     nSA = str_nelem(is->anneal, MAXPTR, ptr1);
3278     if (nSA == 1 && (ptr1[0][0] == 'n' || ptr1[0][0] == 'N'))
3279     {
3280         nSA = 0;
3281     }
3282     if (nSA > 0 && nSA != nr)
3283     {
3284         gmx_fatal(FARGS, "Not enough annealing values: %d (for %d groups)\n", nSA, nr);
3285     }
3286     else
3287     {
3288         snew(ir->opts.annealing, nr);
3289         snew(ir->opts.anneal_npoints, nr);
3290         snew(ir->opts.anneal_time, nr);
3291         snew(ir->opts.anneal_temp, nr);
3292         for (i = 0; i < nr; i++)
3293         {
3294             ir->opts.annealing[i]      = eannNO;
3295             ir->opts.anneal_npoints[i] = 0;
3296             ir->opts.anneal_time[i]    = NULL;
3297             ir->opts.anneal_temp[i]    = NULL;
3298         }
3299         if (nSA > 0)
3300         {
3301             bAnneal = FALSE;
3302             for (i = 0; i < nr; i++)
3303             {
3304                 if (ptr1[i][0] == 'n' || ptr1[i][0] == 'N')
3305                 {
3306                     ir->opts.annealing[i] = eannNO;
3307                 }
3308                 else if (ptr1[i][0] == 's' || ptr1[i][0] == 'S')
3309                 {
3310                     ir->opts.annealing[i] = eannSINGLE;
3311                     bAnneal               = TRUE;
3312                 }
3313                 else if (ptr1[i][0] == 'p' || ptr1[i][0] == 'P')
3314                 {
3315                     ir->opts.annealing[i] = eannPERIODIC;
3316                     bAnneal               = TRUE;
3317                 }
3318             }
3319             if (bAnneal)
3320             {
3321                 /* Read the other fields too */
3322                 nSA_points = str_nelem(is->anneal_npoints, MAXPTR, ptr1);
3323                 if (nSA_points != nSA)
3324                 {
3325                     gmx_fatal(FARGS, "Found %d annealing-npoints values for %d groups\n", nSA_points, nSA);
3326                 }
3327                 for (k = 0, i = 0; i < nr; i++)
3328                 {
3329                     ir->opts.anneal_npoints[i] = strtol(ptr1[i], NULL, 10);
3330                     if (ir->opts.anneal_npoints[i] == 1)
3331                     {
3332                         gmx_fatal(FARGS, "Please specify at least a start and an end point for annealing\n");
3333                     }
3334                     snew(ir->opts.anneal_time[i], ir->opts.anneal_npoints[i]);
3335                     snew(ir->opts.anneal_temp[i], ir->opts.anneal_npoints[i]);
3336                     k += ir->opts.anneal_npoints[i];
3337                 }
3338
3339                 nSA_time = str_nelem(is->anneal_time, MAXPTR, ptr1);
3340                 if (nSA_time != k)
3341                 {
3342                     gmx_fatal(FARGS, "Found %d annealing-time values, wanter %d\n", nSA_time, k);
3343                 }
3344                 nSA_temp = str_nelem(is->anneal_temp, MAXPTR, ptr2);
3345                 if (nSA_temp != k)
3346                 {
3347                     gmx_fatal(FARGS, "Found %d annealing-temp values, wanted %d\n", nSA_temp, k);
3348                 }
3349
3350                 for (i = 0, k = 0; i < nr; i++)
3351                 {
3352
3353                     for (j = 0; j < ir->opts.anneal_npoints[i]; j++)
3354                     {
3355                         ir->opts.anneal_time[i][j] = strtod(ptr1[k], NULL);
3356                         ir->opts.anneal_temp[i][j] = strtod(ptr2[k], NULL);
3357                         if (j == 0)
3358                         {
3359                             if (ir->opts.anneal_time[i][0] > (ir->init_t+GMX_REAL_EPS))
3360                             {
3361                                 gmx_fatal(FARGS, "First time point for annealing > init_t.\n");
3362                             }
3363                         }
3364                         else
3365                         {
3366                             /* j>0 */
3367                             if (ir->opts.anneal_time[i][j] < ir->opts.anneal_time[i][j-1])
3368                             {
3369                                 gmx_fatal(FARGS, "Annealing timepoints out of order: t=%f comes after t=%f\n",
3370                                           ir->opts.anneal_time[i][j], ir->opts.anneal_time[i][j-1]);
3371                             }
3372                         }
3373                         if (ir->opts.anneal_temp[i][j] < 0)
3374                         {
3375                             gmx_fatal(FARGS, "Found negative temperature in annealing: %f\n", ir->opts.anneal_temp[i][j]);
3376                         }
3377                         k++;
3378                     }
3379                 }
3380                 /* Print out some summary information, to make sure we got it right */
3381                 for (i = 0, k = 0; i < nr; i++)
3382                 {
3383                     if (ir->opts.annealing[i] != eannNO)
3384                     {
3385                         j = groups->grps[egcTC].nm_ind[i];
3386                         fprintf(stderr, "Simulated annealing for group %s: %s, %d timepoints\n",
3387                                 *(groups->grpname[j]), eann_names[ir->opts.annealing[i]],
3388                                 ir->opts.anneal_npoints[i]);
3389                         fprintf(stderr, "Time (ps)   Temperature (K)\n");
3390                         /* All terms except the last one */
3391                         for (j = 0; j < (ir->opts.anneal_npoints[i]-1); j++)
3392                         {
3393                             fprintf(stderr, "%9.1f      %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3394                         }
3395
3396                         /* Finally the last one */
3397                         j = ir->opts.anneal_npoints[i]-1;
3398                         if (ir->opts.annealing[i] == eannSINGLE)
3399                         {
3400                             fprintf(stderr, "%9.1f-     %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3401                         }
3402                         else
3403                         {
3404                             fprintf(stderr, "%9.1f      %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3405                             if (fabs(ir->opts.anneal_temp[i][j]-ir->opts.anneal_temp[i][0]) > GMX_REAL_EPS)
3406                             {
3407                                 warning_note(wi, "There is a temperature jump when your annealing loops back.\n");
3408                             }
3409                         }
3410                     }
3411                 }
3412             }
3413         }
3414     }
3415
3416     if (ir->ePull != epullNO)
3417     {
3418         make_pull_groups(ir->pull, is->pull_grp, grps, gnames);
3419
3420         make_pull_coords(ir->pull);
3421     }
3422
3423     if (ir->bRot)
3424     {
3425         make_rotation_groups(ir->rot, is->rot_grp, grps, gnames);
3426     }
3427
3428     if (ir->eSwapCoords != eswapNO)
3429     {
3430         make_swap_groups(ir->swap, swapgrp, splitgrp0, splitgrp1, solgrp, grps, gnames);
3431     }
3432
3433     /* Make indices for IMD session */
3434     if (ir->bIMD)
3435     {
3436         make_IMD_group(ir->imd, is->imd_grp, grps, gnames);
3437     }
3438
3439     nacc = str_nelem(is->acc, MAXPTR, ptr1);
3440     nacg = str_nelem(is->accgrps, MAXPTR, ptr2);
3441     if (nacg*DIM != nacc)
3442     {
3443         gmx_fatal(FARGS, "Invalid Acceleration input: %d groups and %d acc. values",
3444                   nacg, nacc);
3445     }
3446     do_numbering(natoms, groups, nacg, ptr2, grps, gnames, egcACC,
3447                  restnm, egrptpALL_GENREST, bVerbose, wi);
3448     nr = groups->grps[egcACC].nr;
3449     snew(ir->opts.acc, nr);
3450     ir->opts.ngacc = nr;
3451
3452     for (i = k = 0; (i < nacg); i++)
3453     {
3454         for (j = 0; (j < DIM); j++, k++)
3455         {
3456             ir->opts.acc[i][j] = strtod(ptr1[k], NULL);
3457         }
3458     }
3459     for (; (i < nr); i++)
3460     {
3461         for (j = 0; (j < DIM); j++)
3462         {
3463             ir->opts.acc[i][j] = 0;
3464         }
3465     }
3466
3467     nfrdim  = str_nelem(is->frdim, MAXPTR, ptr1);
3468     nfreeze = str_nelem(is->freeze, MAXPTR, ptr2);
3469     if (nfrdim != DIM*nfreeze)
3470     {
3471         gmx_fatal(FARGS, "Invalid Freezing input: %d groups and %d freeze values",
3472                   nfreeze, nfrdim);
3473     }
3474     do_numbering(natoms, groups, nfreeze, ptr2, grps, gnames, egcFREEZE,
3475                  restnm, egrptpALL_GENREST, bVerbose, wi);
3476     nr             = groups->grps[egcFREEZE].nr;
3477     ir->opts.ngfrz = nr;
3478     snew(ir->opts.nFreeze, nr);
3479     for (i = k = 0; (i < nfreeze); i++)
3480     {
3481         for (j = 0; (j < DIM); j++, k++)
3482         {
3483             ir->opts.nFreeze[i][j] = (gmx_strncasecmp(ptr1[k], "Y", 1) == 0);
3484             if (!ir->opts.nFreeze[i][j])
3485             {
3486                 if (gmx_strncasecmp(ptr1[k], "N", 1) != 0)
3487                 {
3488                     sprintf(warnbuf, "Please use Y(ES) or N(O) for freezedim only "
3489                             "(not %s)", ptr1[k]);
3490                     warning(wi, warn_buf);
3491                 }
3492             }
3493         }
3494     }
3495     for (; (i < nr); i++)
3496     {
3497         for (j = 0; (j < DIM); j++)
3498         {
3499             ir->opts.nFreeze[i][j] = 0;
3500         }
3501     }
3502
3503     nenergy = str_nelem(is->energy, MAXPTR, ptr1);
3504     do_numbering(natoms, groups, nenergy, ptr1, grps, gnames, egcENER,
3505                  restnm, egrptpALL_GENREST, bVerbose, wi);
3506     add_wall_energrps(groups, ir->nwall, symtab);
3507     ir->opts.ngener = groups->grps[egcENER].nr;
3508     nvcm            = str_nelem(is->vcm, MAXPTR, ptr1);
3509     bRest           =
3510         do_numbering(natoms, groups, nvcm, ptr1, grps, gnames, egcVCM,
3511                      restnm, nvcm == 0 ? egrptpALL_GENREST : egrptpPART, bVerbose, wi);
3512     if (bRest)
3513     {
3514         warning(wi, "Some atoms are not part of any center of mass motion removal group.\n"
3515                 "This may lead to artifacts.\n"
3516                 "In most cases one should use one group for the whole system.");
3517     }
3518
3519     /* Now we have filled the freeze struct, so we can calculate NRDF */
3520     calc_nrdf(mtop, ir, gnames);
3521
3522     if (v && NULL)
3523     {
3524         real fac, ntot = 0;
3525
3526         /* Must check per group! */
3527         for (i = 0; (i < ir->opts.ngtc); i++)
3528         {
3529             ntot += ir->opts.nrdf[i];
3530         }
3531         if (ntot != (DIM*natoms))
3532         {
3533             fac = sqrt(ntot/(DIM*natoms));
3534             if (bVerbose)
3535             {
3536                 fprintf(stderr, "Scaling velocities by a factor of %.3f to account for constraints\n"
3537                         "and removal of center of mass motion\n", fac);
3538             }
3539             for (i = 0; (i < natoms); i++)
3540             {
3541                 svmul(fac, v[i], v[i]);
3542             }
3543         }
3544     }
3545
3546     nuser = str_nelem(is->user1, MAXPTR, ptr1);
3547     do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcUser1,
3548                  restnm, egrptpALL_GENREST, bVerbose, wi);
3549     nuser = str_nelem(is->user2, MAXPTR, ptr1);
3550     do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcUser2,
3551                  restnm, egrptpALL_GENREST, bVerbose, wi);
3552     nuser = str_nelem(is->x_compressed_groups, MAXPTR, ptr1);
3553     do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcCompressedX,
3554                  restnm, egrptpONE, bVerbose, wi);
3555     nofg = str_nelem(is->orirefitgrp, MAXPTR, ptr1);
3556     do_numbering(natoms, groups, nofg, ptr1, grps, gnames, egcORFIT,
3557                  restnm, egrptpALL_GENREST, bVerbose, wi);
3558
3559     /* QMMM input processing */
3560     nQMg          = str_nelem(is->QMMM, MAXPTR, ptr1);
3561     nQMmethod     = str_nelem(is->QMmethod, MAXPTR, ptr2);
3562     nQMbasis      = str_nelem(is->QMbasis, MAXPTR, ptr3);
3563     if ((nQMmethod != nQMg) || (nQMbasis != nQMg))
3564     {
3565         gmx_fatal(FARGS, "Invalid QMMM input: %d groups %d basissets"
3566                   " and %d methods\n", nQMg, nQMbasis, nQMmethod);
3567     }
3568     /* group rest, if any, is always MM! */
3569     do_numbering(natoms, groups, nQMg, ptr1, grps, gnames, egcQMMM,
3570                  restnm, egrptpALL_GENREST, bVerbose, wi);
3571     nr            = nQMg; /*atoms->grps[egcQMMM].nr;*/
3572     ir->opts.ngQM = nQMg;
3573     snew(ir->opts.QMmethod, nr);
3574     snew(ir->opts.QMbasis, nr);
3575     for (i = 0; i < nr; i++)
3576     {
3577         /* input consists of strings: RHF CASSCF PM3 .. These need to be
3578          * converted to the corresponding enum in names.c
3579          */
3580         ir->opts.QMmethod[i] = search_QMstring(ptr2[i], eQMmethodNR,
3581                                                eQMmethod_names);
3582         ir->opts.QMbasis[i]  = search_QMstring(ptr3[i], eQMbasisNR,
3583                                                eQMbasis_names);
3584
3585     }
3586     nQMmult   = str_nelem(is->QMmult, MAXPTR, ptr1);
3587     nQMcharge = str_nelem(is->QMcharge, MAXPTR, ptr2);
3588     nbSH      = str_nelem(is->bSH, MAXPTR, ptr3);
3589     snew(ir->opts.QMmult, nr);
3590     snew(ir->opts.QMcharge, nr);
3591     snew(ir->opts.bSH, nr);
3592
3593     for (i = 0; i < nr; i++)
3594     {
3595         ir->opts.QMmult[i]   = strtol(ptr1[i], NULL, 10);
3596         ir->opts.QMcharge[i] = strtol(ptr2[i], NULL, 10);
3597         ir->opts.bSH[i]      = (gmx_strncasecmp(ptr3[i], "Y", 1) == 0);
3598     }
3599
3600     nCASelec  = str_nelem(is->CASelectrons, MAXPTR, ptr1);
3601     nCASorb   = str_nelem(is->CASorbitals, MAXPTR, ptr2);
3602     snew(ir->opts.CASelectrons, nr);
3603     snew(ir->opts.CASorbitals, nr);
3604     for (i = 0; i < nr; i++)
3605     {
3606         ir->opts.CASelectrons[i] = strtol(ptr1[i], NULL, 10);
3607         ir->opts.CASorbitals[i]  = strtol(ptr2[i], NULL, 10);
3608     }
3609     /* special optimization options */
3610
3611     nbOPT = str_nelem(is->bOPT, MAXPTR, ptr1);
3612     nbTS  = str_nelem(is->bTS, MAXPTR, ptr2);
3613     snew(ir->opts.bOPT, nr);
3614     snew(ir->opts.bTS, nr);
3615     for (i = 0; i < nr; i++)
3616     {
3617         ir->opts.bOPT[i] = (gmx_strncasecmp(ptr1[i], "Y", 1) == 0);
3618         ir->opts.bTS[i]  = (gmx_strncasecmp(ptr2[i], "Y", 1) == 0);
3619     }
3620     nSAon     = str_nelem(is->SAon, MAXPTR, ptr1);
3621     nSAoff    = str_nelem(is->SAoff, MAXPTR, ptr2);
3622     nSAsteps  = str_nelem(is->SAsteps, MAXPTR, ptr3);
3623     snew(ir->opts.SAon, nr);
3624     snew(ir->opts.SAoff, nr);
3625     snew(ir->opts.SAsteps, nr);
3626
3627     for (i = 0; i < nr; i++)
3628     {
3629         ir->opts.SAon[i]    = strtod(ptr1[i], NULL);
3630         ir->opts.SAoff[i]   = strtod(ptr2[i], NULL);
3631         ir->opts.SAsteps[i] = strtol(ptr3[i], NULL, 10);
3632     }
3633     /* end of QMMM input */
3634
3635     if (bVerbose)
3636     {
3637         for (i = 0; (i < egcNR); i++)
3638         {
3639             fprintf(stderr, "%-16s has %d element(s):", gtypes[i], groups->grps[i].nr);
3640             for (j = 0; (j < groups->grps[i].nr); j++)
3641             {
3642                 fprintf(stderr, " %s", *(groups->grpname[groups->grps[i].nm_ind[j]]));
3643             }
3644             fprintf(stderr, "\n");
3645         }
3646     }
3647
3648     nr = groups->grps[egcENER].nr;
3649     snew(ir->opts.egp_flags, nr*nr);
3650
3651     bExcl = do_egp_flag(ir, groups, "energygrp-excl", is->egpexcl, EGP_EXCL);
3652     if (bExcl && ir->cutoff_scheme == ecutsVERLET)
3653     {
3654         warning_error(wi, "Energy group exclusions are not (yet) implemented for the Verlet scheme");
3655     }
3656     if (bExcl && EEL_FULL(ir->coulombtype))
3657     {
3658         warning(wi, "Can not exclude the lattice Coulomb energy between energy groups");
3659     }
3660
3661     bTable = do_egp_flag(ir, groups, "energygrp-table", is->egptable, EGP_TABLE);
3662     if (bTable && !(ir->vdwtype == evdwUSER) &&
3663         !(ir->coulombtype == eelUSER) && !(ir->coulombtype == eelPMEUSER) &&
3664         !(ir->coulombtype == eelPMEUSERSWITCH))
3665     {
3666         gmx_fatal(FARGS, "Can only have energy group pair tables in combination with user tables for VdW and/or Coulomb");
3667     }
3668
3669     decode_cos(is->efield_x, &(ir->ex[XX]));
3670     decode_cos(is->efield_xt, &(ir->et[XX]));
3671     decode_cos(is->efield_y, &(ir->ex[YY]));
3672     decode_cos(is->efield_yt, &(ir->et[YY]));
3673     decode_cos(is->efield_z, &(ir->ex[ZZ]));
3674     decode_cos(is->efield_zt, &(ir->et[ZZ]));
3675
3676     if (ir->bAdress)
3677     {
3678         do_adress_index(ir->adress, groups, gnames, &(ir->opts), wi);
3679     }
3680
3681     for (i = 0; (i < grps->nr); i++)
3682     {
3683         sfree(gnames[i]);
3684     }
3685     sfree(gnames);
3686     done_blocka(grps);
3687     sfree(grps);
3688
3689 }
3690
3691
3692
3693 static void check_disre(gmx_mtop_t *mtop)
3694 {
3695     gmx_ffparams_t *ffparams;
3696     t_functype     *functype;
3697     t_iparams      *ip;
3698     int             i, ndouble, ftype;
3699     int             label, old_label;
3700
3701     if (gmx_mtop_ftype_count(mtop, F_DISRES) > 0)
3702     {
3703         ffparams  = &mtop->ffparams;
3704         functype  = ffparams->functype;
3705         ip        = ffparams->iparams;
3706         ndouble   = 0;
3707         old_label = -1;
3708         for (i = 0; i < ffparams->ntypes; i++)
3709         {
3710             ftype = functype[i];
3711             if (ftype == F_DISRES)
3712             {
3713                 label = ip[i].disres.label;
3714                 if (label == old_label)
3715                 {
3716                     fprintf(stderr, "Distance restraint index %d occurs twice\n", label);
3717                     ndouble++;
3718                 }
3719                 old_label = label;
3720             }
3721         }
3722         if (ndouble > 0)
3723         {
3724             gmx_fatal(FARGS, "Found %d double distance restraint indices,\n"
3725                       "probably the parameters for multiple pairs in one restraint "
3726                       "are not identical\n", ndouble);
3727         }
3728     }
3729 }
3730
3731 static gmx_bool absolute_reference(t_inputrec *ir, gmx_mtop_t *sys,
3732                                    gmx_bool posres_only,
3733                                    ivec AbsRef)
3734 {
3735     int                  d, g, i;
3736     gmx_mtop_ilistloop_t iloop;
3737     t_ilist             *ilist;
3738     int                  nmol;
3739     t_iparams           *pr;
3740
3741     clear_ivec(AbsRef);
3742
3743     if (!posres_only)
3744     {
3745         /* Check the COM */
3746         for (d = 0; d < DIM; d++)
3747         {
3748             AbsRef[d] = (d < ndof_com(ir) ? 0 : 1);
3749         }
3750         /* Check for freeze groups */
3751         for (g = 0; g < ir->opts.ngfrz; g++)
3752         {
3753             for (d = 0; d < DIM; d++)
3754             {
3755                 if (ir->opts.nFreeze[g][d] != 0)
3756                 {
3757                     AbsRef[d] = 1;
3758                 }
3759             }
3760         }
3761     }
3762
3763     /* Check for position restraints */
3764     iloop = gmx_mtop_ilistloop_init(sys);
3765     while (gmx_mtop_ilistloop_next(iloop, &ilist, &nmol))
3766     {
3767         if (nmol > 0 &&
3768             (AbsRef[XX] == 0 || AbsRef[YY] == 0 || AbsRef[ZZ] == 0))
3769         {
3770             for (i = 0; i < ilist[F_POSRES].nr; i += 2)
3771             {
3772                 pr = &sys->ffparams.iparams[ilist[F_POSRES].iatoms[i]];
3773                 for (d = 0; d < DIM; d++)
3774                 {
3775                     if (pr->posres.fcA[d] != 0)
3776                     {
3777                         AbsRef[d] = 1;
3778                     }
3779                 }
3780             }
3781             for (i = 0; i < ilist[F_FBPOSRES].nr; i += 2)
3782             {
3783                 /* Check for flat-bottom posres */
3784                 pr = &sys->ffparams.iparams[ilist[F_FBPOSRES].iatoms[i]];
3785                 if (pr->fbposres.k != 0)
3786                 {
3787                     switch (pr->fbposres.geom)
3788                     {
3789                         case efbposresSPHERE:
3790                             AbsRef[XX] = AbsRef[YY] = AbsRef[ZZ] = 1;
3791                             break;
3792                         case efbposresCYLINDER:
3793                             AbsRef[XX] = AbsRef[YY] = 1;
3794                             break;
3795                         case efbposresX: /* d=XX */
3796                         case efbposresY: /* d=YY */
3797                         case efbposresZ: /* d=ZZ */
3798                             d         = pr->fbposres.geom - efbposresX;
3799                             AbsRef[d] = 1;
3800                             break;
3801                         default:
3802                             gmx_fatal(FARGS, " Invalid geometry for flat-bottom position restraint.\n"
3803                                       "Expected nr between 1 and %d. Found %d\n", efbposresNR-1,
3804                                       pr->fbposres.geom);
3805                     }
3806                 }
3807             }
3808         }
3809     }
3810
3811     return (AbsRef[XX] != 0 && AbsRef[YY] != 0 && AbsRef[ZZ] != 0);
3812 }
3813
3814 static void
3815 check_combination_rule_differences(const gmx_mtop_t *mtop, int state,
3816                                    gmx_bool *bC6ParametersWorkWithGeometricRules,
3817                                    gmx_bool *bC6ParametersWorkWithLBRules,
3818                                    gmx_bool *bLBRulesPossible)
3819 {
3820     int           ntypes, tpi, tpj, thisLBdiff, thisgeomdiff;
3821     int          *typecount;
3822     real          tol;
3823     double        geometricdiff, LBdiff;
3824     double        c6i, c6j, c12i, c12j;
3825     double        c6, c6_geometric, c6_LB;
3826     double        sigmai, sigmaj, epsi, epsj;
3827     gmx_bool      bCanDoLBRules, bCanDoGeometricRules;
3828     const char   *ptr;
3829
3830     /* A tolerance of 1e-5 seems reasonable for (possibly hand-typed)
3831      * force-field floating point parameters.
3832      */
3833     tol = 1e-5;
3834     ptr = getenv("GMX_LJCOMB_TOL");
3835     if (ptr != NULL)
3836     {
3837         double dbl;
3838
3839         sscanf(ptr, "%lf", &dbl);
3840         tol = dbl;
3841     }
3842
3843     *bC6ParametersWorkWithLBRules         = TRUE;
3844     *bC6ParametersWorkWithGeometricRules  = TRUE;
3845     bCanDoLBRules                         = TRUE;
3846     bCanDoGeometricRules                  = TRUE;
3847     ntypes                                = mtop->ffparams.atnr;
3848     snew(typecount, ntypes);
3849     gmx_mtop_count_atomtypes(mtop, state, typecount);
3850     geometricdiff           = LBdiff = 0.0;
3851     *bLBRulesPossible       = TRUE;
3852     for (tpi = 0; tpi < ntypes; ++tpi)
3853     {
3854         c6i  = mtop->ffparams.iparams[(ntypes + 1) * tpi].lj.c6;
3855         c12i = mtop->ffparams.iparams[(ntypes + 1) * tpi].lj.c12;
3856         for (tpj = tpi; tpj < ntypes; ++tpj)
3857         {
3858             c6j          = mtop->ffparams.iparams[(ntypes + 1) * tpj].lj.c6;
3859             c12j         = mtop->ffparams.iparams[(ntypes + 1) * tpj].lj.c12;
3860             c6           = mtop->ffparams.iparams[ntypes * tpi + tpj].lj.c6;
3861             c6_geometric = sqrt(c6i * c6j);
3862             if (!gmx_numzero(c6_geometric))
3863             {
3864                 if (!gmx_numzero(c12i) && !gmx_numzero(c12j))
3865                 {
3866                     sigmai   = pow(c12i / c6i, 1.0/6.0);
3867                     sigmaj   = pow(c12j / c6j, 1.0/6.0);
3868                     epsi     = c6i * c6i /(4.0 * c12i);
3869                     epsj     = c6j * c6j /(4.0 * c12j);
3870                     c6_LB    = 4.0 * pow(epsi * epsj, 1.0/2.0) * pow(0.5 * (sigmai + sigmaj), 6);
3871                 }
3872                 else
3873                 {
3874                     *bLBRulesPossible = FALSE;
3875                     c6_LB             = c6_geometric;
3876                 }
3877                 bCanDoLBRules = gmx_within_tol(c6_LB, c6, tol);
3878             }
3879
3880             if (FALSE == bCanDoLBRules)
3881             {
3882                 *bC6ParametersWorkWithLBRules = FALSE;
3883             }
3884
3885             bCanDoGeometricRules = gmx_within_tol(c6_geometric, c6, tol);
3886
3887             if (FALSE == bCanDoGeometricRules)
3888             {
3889                 *bC6ParametersWorkWithGeometricRules = FALSE;
3890             }
3891         }
3892     }
3893     sfree(typecount);
3894 }
3895
3896 static void
3897 check_combination_rules(const t_inputrec *ir, const gmx_mtop_t *mtop,
3898                         warninp_t wi)
3899 {
3900     char     err_buf[256];
3901     gmx_bool bLBRulesPossible, bC6ParametersWorkWithGeometricRules, bC6ParametersWorkWithLBRules;
3902
3903     check_combination_rule_differences(mtop, 0,
3904                                        &bC6ParametersWorkWithGeometricRules,
3905                                        &bC6ParametersWorkWithLBRules,
3906                                        &bLBRulesPossible);
3907     if (ir->ljpme_combination_rule == eljpmeLB)
3908     {
3909         if (FALSE == bC6ParametersWorkWithLBRules || FALSE == bLBRulesPossible)
3910         {
3911             warning(wi, "You are using arithmetic-geometric combination rules "
3912                     "in LJ-PME, but your non-bonded C6 parameters do not "
3913                     "follow these rules.");
3914         }
3915     }
3916     else
3917     {
3918         if (FALSE == bC6ParametersWorkWithGeometricRules)
3919         {
3920             if (ir->eDispCorr != edispcNO)
3921             {
3922                 warning_note(wi, "You are using geometric combination rules in "
3923                              "LJ-PME, but your non-bonded C6 parameters do "
3924                              "not follow these rules. "
3925                              "This will introduce very small errors in the forces and energies in "
3926                              "your simulations. Dispersion correction will correct total energy "
3927                              "and/or pressure for isotropic systems, but not forces or surface tensions.");
3928             }
3929             else
3930             {
3931                 warning_note(wi, "You are using geometric combination rules in "
3932                              "LJ-PME, but your non-bonded C6 parameters do "
3933                              "not follow these rules. "
3934                              "This will introduce very small errors in the forces and energies in "
3935                              "your simulations. If your system is homogeneous, consider using dispersion correction "
3936                              "for the total energy and pressure.");
3937             }
3938         }
3939     }
3940 }
3941
3942 void triple_check(const char *mdparin, t_inputrec *ir, gmx_mtop_t *sys,
3943                   warninp_t wi)
3944 {
3945     char                      err_buf[STRLEN];
3946     int                       i, m, c, nmol, npct;
3947     gmx_bool                  bCharge, bAcc;
3948     real                      gdt_max, *mgrp, mt;
3949     rvec                      acc;
3950     gmx_mtop_atomloop_block_t aloopb;
3951     gmx_mtop_atomloop_all_t   aloop;
3952     t_atom                   *atom;
3953     ivec                      AbsRef;
3954     char                      warn_buf[STRLEN];
3955
3956     set_warning_line(wi, mdparin, -1);
3957
3958     if (ir->cutoff_scheme == ecutsVERLET &&
3959         ir->verletbuf_tol > 0 &&
3960         ir->nstlist > 1 &&
3961         ((EI_MD(ir->eI) || EI_SD(ir->eI)) &&
3962          (ir->etc == etcVRESCALE || ir->etc == etcBERENDSEN)))
3963     {
3964         /* Check if a too small Verlet buffer might potentially
3965          * cause more drift than the thermostat can couple off.
3966          */
3967         /* Temperature error fraction for warning and suggestion */
3968         const real T_error_warn    = 0.002;
3969         const real T_error_suggest = 0.001;
3970         /* For safety: 2 DOF per atom (typical with constraints) */
3971         const real nrdf_at         = 2;
3972         real       T, tau, max_T_error;
3973         int        i;
3974
3975         T   = 0;
3976         tau = 0;
3977         for (i = 0; i < ir->opts.ngtc; i++)
3978         {
3979             T   = max(T, ir->opts.ref_t[i]);
3980             tau = max(tau, ir->opts.tau_t[i]);
3981         }
3982         if (T > 0)
3983         {
3984             /* This is a worst case estimate of the temperature error,
3985              * assuming perfect buffer estimation and no cancelation
3986              * of errors. The factor 0.5 is because energy distributes
3987              * equally over Ekin and Epot.
3988              */
3989             max_T_error = 0.5*tau*ir->verletbuf_tol/(nrdf_at*BOLTZ*T);
3990             if (max_T_error > T_error_warn)
3991             {
3992                 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.",
3993                         ir->verletbuf_tol, T, tau,
3994                         100*max_T_error,
3995                         100*T_error_suggest,
3996                         ir->verletbuf_tol*T_error_suggest/max_T_error);
3997                 warning(wi, warn_buf);
3998             }
3999         }
4000     }
4001
4002     if (ETC_ANDERSEN(ir->etc))
4003     {
4004         int i;
4005
4006         for (i = 0; i < ir->opts.ngtc; i++)
4007         {
4008             sprintf(err_buf, "all tau_t must currently be equal using Andersen temperature control, violated for group %d", i);
4009             CHECK(ir->opts.tau_t[0] != ir->opts.tau_t[i]);
4010             sprintf(err_buf, "all tau_t must be postive using Andersen temperature control, tau_t[%d]=%10.6f",
4011                     i, ir->opts.tau_t[i]);
4012             CHECK(ir->opts.tau_t[i] < 0);
4013         }
4014
4015         for (i = 0; i < ir->opts.ngtc; i++)
4016         {
4017             int nsteps = (int)(ir->opts.tau_t[i]/ir->delta_t);
4018             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);
4019             CHECK((nsteps % ir->nstcomm) && (ir->etc == etcANDERSENMASSIVE));
4020         }
4021     }
4022
4023     if (EI_DYNAMICS(ir->eI) && !EI_SD(ir->eI) && ir->eI != eiBD &&
4024         ir->comm_mode == ecmNO &&
4025         !(absolute_reference(ir, sys, FALSE, AbsRef) || ir->nsteps <= 10) &&
4026         !ETC_ANDERSEN(ir->etc))
4027     {
4028         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");
4029     }
4030
4031     /* Check for pressure coupling with absolute position restraints */
4032     if (ir->epc != epcNO && ir->refcoord_scaling == erscNO)
4033     {
4034         absolute_reference(ir, sys, TRUE, AbsRef);
4035         {
4036             for (m = 0; m < DIM; m++)
4037             {
4038                 if (AbsRef[m] && norm2(ir->compress[m]) > 0)
4039                 {
4040                     warning(wi, "You are using pressure coupling with absolute position restraints, this will give artifacts. Use the refcoord_scaling option.");
4041                     break;
4042                 }
4043             }
4044         }
4045     }
4046
4047     bCharge = FALSE;
4048     aloopb  = gmx_mtop_atomloop_block_init(sys);
4049     while (gmx_mtop_atomloop_block_next(aloopb, &atom, &nmol))
4050     {
4051         if (atom->q != 0 || atom->qB != 0)
4052         {
4053             bCharge = TRUE;
4054         }
4055     }
4056
4057     if (!bCharge)
4058     {
4059         if (EEL_FULL(ir->coulombtype))
4060         {
4061             sprintf(err_buf,
4062                     "You are using full electrostatics treatment %s for a system without charges.\n"
4063                     "This costs a lot of performance for just processing zeros, consider using %s instead.\n",
4064                     EELTYPE(ir->coulombtype), EELTYPE(eelCUT));
4065             warning(wi, err_buf);
4066         }
4067     }
4068     else
4069     {
4070         if (ir->coulombtype == eelCUT && ir->rcoulomb > 0 && !ir->implicit_solvent)
4071         {
4072             sprintf(err_buf,
4073                     "You are using a plain Coulomb cut-off, which might produce artifacts.\n"
4074                     "You might want to consider using %s electrostatics.\n",
4075                     EELTYPE(eelPME));
4076             warning_note(wi, err_buf);
4077         }
4078     }
4079
4080     /* Check if combination rules used in LJ-PME are the same as in the force field */
4081     if (EVDW_PME(ir->vdwtype))
4082     {
4083         check_combination_rules(ir, sys, wi);
4084     }
4085
4086     /* Generalized reaction field */
4087     if (ir->opts.ngtc == 0)
4088     {
4089         sprintf(err_buf, "No temperature coupling while using coulombtype %s",
4090                 eel_names[eelGRF]);
4091         CHECK(ir->coulombtype == eelGRF);
4092     }
4093     else
4094     {
4095         sprintf(err_buf, "When using coulombtype = %s"
4096                 " ref-t for temperature coupling should be > 0",
4097                 eel_names[eelGRF]);
4098         CHECK((ir->coulombtype == eelGRF) && (ir->opts.ref_t[0] <= 0));
4099     }
4100
4101     if (ir->eI == eiSD1 &&
4102         (gmx_mtop_ftype_count(sys, F_CONSTR) > 0 ||
4103          gmx_mtop_ftype_count(sys, F_SETTLE) > 0))
4104     {
4105         sprintf(warn_buf, "With constraints integrator %s is less accurate, consider using %s instead", ei_names[ir->eI], ei_names[eiSD2]);
4106         warning_note(wi, warn_buf);
4107     }
4108
4109     bAcc = FALSE;
4110     for (i = 0; (i < sys->groups.grps[egcACC].nr); i++)
4111     {
4112         for (m = 0; (m < DIM); m++)
4113         {
4114             if (fabs(ir->opts.acc[i][m]) > 1e-6)
4115             {
4116                 bAcc = TRUE;
4117             }
4118         }
4119     }
4120     if (bAcc)
4121     {
4122         clear_rvec(acc);
4123         snew(mgrp, sys->groups.grps[egcACC].nr);
4124         aloop = gmx_mtop_atomloop_all_init(sys);
4125         while (gmx_mtop_atomloop_all_next(aloop, &i, &atom))
4126         {
4127             mgrp[ggrpnr(&sys->groups, egcACC, i)] += atom->m;
4128         }
4129         mt = 0.0;
4130         for (i = 0; (i < sys->groups.grps[egcACC].nr); i++)
4131         {
4132             for (m = 0; (m < DIM); m++)
4133             {
4134                 acc[m] += ir->opts.acc[i][m]*mgrp[i];
4135             }
4136             mt += mgrp[i];
4137         }
4138         for (m = 0; (m < DIM); m++)
4139         {
4140             if (fabs(acc[m]) > 1e-6)
4141             {
4142                 const char *dim[DIM] = { "X", "Y", "Z" };
4143                 fprintf(stderr,
4144                         "Net Acceleration in %s direction, will %s be corrected\n",
4145                         dim[m], ir->nstcomm != 0 ? "" : "not");
4146                 if (ir->nstcomm != 0 && m < ndof_com(ir))
4147                 {
4148                     acc[m] /= mt;
4149                     for (i = 0; (i < sys->groups.grps[egcACC].nr); i++)
4150                     {
4151                         ir->opts.acc[i][m] -= acc[m];
4152                     }
4153                 }
4154             }
4155         }
4156         sfree(mgrp);
4157     }
4158
4159     if (ir->efep != efepNO && ir->fepvals->sc_alpha != 0 &&
4160         !gmx_within_tol(sys->ffparams.reppow, 12.0, 10*GMX_DOUBLE_EPS))
4161     {
4162         gmx_fatal(FARGS, "Soft-core interactions are only supported with VdW repulsion power 12");
4163     }
4164
4165     if (ir->ePull != epullNO)
4166     {
4167         gmx_bool bPullAbsoluteRef;
4168
4169         bPullAbsoluteRef = FALSE;
4170         for (i = 0; i < ir->pull->ncoord; i++)
4171         {
4172             bPullAbsoluteRef = bPullAbsoluteRef ||
4173                 ir->pull->coord[i].group[0] == 0 ||
4174                 ir->pull->coord[i].group[1] == 0;
4175         }
4176         if (bPullAbsoluteRef)
4177         {
4178             absolute_reference(ir, sys, FALSE, AbsRef);
4179             for (m = 0; m < DIM; m++)
4180             {
4181                 if (ir->pull->dim[m] && !AbsRef[m])
4182                 {
4183                     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.");
4184                     break;
4185                 }
4186             }
4187         }
4188
4189         if (ir->pull->eGeom == epullgDIRPBC)
4190         {
4191             for (i = 0; i < 3; i++)
4192             {
4193                 for (m = 0; m <= i; m++)
4194                 {
4195                     if ((ir->epc != epcNO && ir->compress[i][m] != 0) ||
4196                         ir->deform[i][m] != 0)
4197                     {
4198                         for (c = 0; c < ir->pull->ncoord; c++)
4199                         {
4200                             if (ir->pull->coord[c].vec[m] != 0)
4201                             {
4202                                 gmx_fatal(FARGS, "Can not have dynamic box while using pull geometry '%s' (dim %c)", EPULLGEOM(ir->pull->eGeom), 'x'+m);
4203                             }
4204                         }
4205                     }
4206                 }
4207             }
4208         }
4209     }
4210
4211     check_disre(sys);
4212 }
4213
4214 void double_check(t_inputrec *ir, matrix box, gmx_bool bConstr, warninp_t wi)
4215 {
4216     real        min_size;
4217     gmx_bool    bTWIN;
4218     char        warn_buf[STRLEN];
4219     const char *ptr;
4220
4221     ptr = check_box(ir->ePBC, box);
4222     if (ptr)
4223     {
4224         warning_error(wi, ptr);
4225     }
4226
4227     if (bConstr && ir->eConstrAlg == econtSHAKE)
4228     {
4229         if (ir->shake_tol <= 0.0)
4230         {
4231             sprintf(warn_buf, "ERROR: shake-tol must be > 0 instead of %g\n",
4232                     ir->shake_tol);
4233             warning_error(wi, warn_buf);
4234         }
4235
4236         if (IR_TWINRANGE(*ir) && ir->nstlist > 1)
4237         {
4238             sprintf(warn_buf, "With twin-range cut-off's and SHAKE the virial and the pressure are incorrect.");
4239             if (ir->epc == epcNO)
4240             {
4241                 warning(wi, warn_buf);
4242             }
4243             else
4244             {
4245                 warning_error(wi, warn_buf);
4246             }
4247         }
4248     }
4249
4250     if ( (ir->eConstrAlg == econtLINCS) && bConstr)
4251     {
4252         /* If we have Lincs constraints: */
4253         if (ir->eI == eiMD && ir->etc == etcNO &&
4254             ir->eConstrAlg == econtLINCS && ir->nLincsIter == 1)
4255         {
4256             sprintf(warn_buf, "For energy conservation with LINCS, lincs_iter should be 2 or larger.\n");
4257             warning_note(wi, warn_buf);
4258         }
4259
4260         if ((ir->eI == eiCG || ir->eI == eiLBFGS) && (ir->nProjOrder < 8))
4261         {
4262             sprintf(warn_buf, "For accurate %s with LINCS constraints, lincs-order should be 8 or more.", ei_names[ir->eI]);
4263             warning_note(wi, warn_buf);
4264         }
4265         if (ir->epc == epcMTTK)
4266         {
4267             warning_error(wi, "MTTK not compatible with lincs -- use shake instead.");
4268         }
4269     }
4270
4271     if (bConstr && ir->epc == epcMTTK)
4272     {
4273         warning_note(wi, "MTTK with constraints is deprecated, and will be removed in GROMACS 5.1");
4274     }
4275
4276     if (ir->LincsWarnAngle > 90.0)
4277     {
4278         sprintf(warn_buf, "lincs-warnangle can not be larger than 90 degrees, setting it to 90.\n");
4279         warning(wi, warn_buf);
4280         ir->LincsWarnAngle = 90.0;
4281     }
4282
4283     if (ir->ePBC != epbcNONE)
4284     {
4285         if (ir->nstlist == 0)
4286         {
4287             warning(wi, "With nstlist=0 atoms are only put into the box at step 0, therefore drifting atoms might cause the simulation to crash.");
4288         }
4289         bTWIN = (ir->rlistlong > ir->rlist);
4290         if (ir->ns_type == ensGRID)
4291         {
4292             if (sqr(ir->rlistlong) >= max_cutoff2(ir->ePBC, box))
4293             {
4294                 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",
4295                         bTWIN ? (ir->rcoulomb == ir->rlistlong ? "rcoulomb" : "rvdw") : "rlist");
4296                 warning_error(wi, warn_buf);
4297             }
4298         }
4299         else
4300         {
4301             min_size = min(box[XX][XX], min(box[YY][YY], box[ZZ][ZZ]));
4302             if (2*ir->rlistlong >= min_size)
4303             {
4304                 sprintf(warn_buf, "ERROR: One of the box lengths is smaller than twice the cut-off length. Increase the box size or decrease rlist.");
4305                 warning_error(wi, warn_buf);
4306                 if (TRICLINIC(box))
4307                 {
4308                     fprintf(stderr, "Grid search might allow larger cut-off's than simple search with triclinic boxes.");
4309                 }
4310             }
4311         }
4312     }
4313 }
4314
4315 void check_chargegroup_radii(const gmx_mtop_t *mtop, const t_inputrec *ir,
4316                              rvec *x,
4317                              warninp_t wi)
4318 {
4319     real rvdw1, rvdw2, rcoul1, rcoul2;
4320     char warn_buf[STRLEN];
4321
4322     calc_chargegroup_radii(mtop, x, &rvdw1, &rvdw2, &rcoul1, &rcoul2);
4323
4324     if (rvdw1 > 0)
4325     {
4326         printf("Largest charge group radii for Van der Waals: %5.3f, %5.3f nm\n",
4327                rvdw1, rvdw2);
4328     }
4329     if (rcoul1 > 0)
4330     {
4331         printf("Largest charge group radii for Coulomb:       %5.3f, %5.3f nm\n",
4332                rcoul1, rcoul2);
4333     }
4334
4335     if (ir->rlist > 0)
4336     {
4337         if (rvdw1  + rvdw2  > ir->rlist ||
4338             rcoul1 + rcoul2 > ir->rlist)
4339         {
4340             sprintf(warn_buf,
4341                     "The sum of the two largest charge group radii (%f) "
4342                     "is larger than rlist (%f)\n",
4343                     max(rvdw1+rvdw2, rcoul1+rcoul2), ir->rlist);
4344             warning(wi, warn_buf);
4345         }
4346         else
4347         {
4348             /* Here we do not use the zero at cut-off macro,
4349              * since user defined interactions might purposely
4350              * not be zero at the cut-off.
4351              */
4352             if (ir_vdw_is_zero_at_cutoff(ir) &&
4353                 rvdw1 + rvdw2 > ir->rlistlong - ir->rvdw)
4354             {
4355                 sprintf(warn_buf, "The sum of the two largest charge group "
4356                         "radii (%f) is larger than %s (%f) - rvdw (%f).\n"
4357                         "With exact cut-offs, better performance can be "
4358                         "obtained with cutoff-scheme = %s, because it "
4359                         "does not use charge groups at all.",
4360                         rvdw1+rvdw2,
4361                         ir->rlistlong > ir->rlist ? "rlistlong" : "rlist",
4362                         ir->rlistlong, ir->rvdw,
4363                         ecutscheme_names[ecutsVERLET]);
4364                 if (ir_NVE(ir))
4365                 {
4366                     warning(wi, warn_buf);
4367                 }
4368                 else
4369                 {
4370                     warning_note(wi, warn_buf);
4371                 }
4372             }
4373             if (ir_coulomb_is_zero_at_cutoff(ir) &&
4374                 rcoul1 + rcoul2 > ir->rlistlong - ir->rcoulomb)
4375             {
4376                 sprintf(warn_buf, "The sum of the two largest charge group radii (%f) is larger than %s (%f) - rcoulomb (%f).\n"
4377                         "With exact cut-offs, better performance can be obtained with cutoff-scheme = %s, because it does not use charge groups at all.",
4378                         rcoul1+rcoul2,
4379                         ir->rlistlong > ir->rlist ? "rlistlong" : "rlist",
4380                         ir->rlistlong, ir->rcoulomb,
4381                         ecutscheme_names[ecutsVERLET]);
4382                 if (ir_NVE(ir))
4383                 {
4384                     warning(wi, warn_buf);
4385                 }
4386                 else
4387                 {
4388                     warning_note(wi, warn_buf);
4389                 }
4390             }
4391         }
4392     }
4393 }