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