src/gromacs/gmxpreprocess/readir.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team.
   6  * Copyright (c) 2013,2014, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 #ifdef HAVE_CONFIG_H
  38 #include <config.h>
  39 #endif
  40
  41 #include <ctype.h>
  42 #include <stdlib.h>
  43 #include <limits.h>
  44 #include "sysstuff.h"
  45 #include "gromacs/utility/smalloc.h"
  46 #include "typedefs.h"
  47 #include "physics.h"
  48 #include "names.h"
  49 #include "gmx_fatal.h"
  50 #include "macros.h"
  51 #include "index.h"
  52 #include "symtab.h"
  53 #include "gromacs/utility/cstringutil.h"
  54 #include "readinp.h"
  55 #include "warninp.h"
  56 #include "readir.h"
  57 #include "toputil.h"
  58 #include "index.h"
  59 #include "network.h"
  60 #include "vec.h"
  61 #include "pbc.h"
  62 #include "mtop_util.h"
  63 #include "chargegroup.h"
  64 #include "inputrec.h"
  65 #include "calc_verletbuf.h"
  66
  67 #define MAXPTR 254
  68 #define NOGID  255
  69
  70 /* Resource parameters
  71  * Do not change any of these until you read the instruction
  72  * in readinp.h. Some cpp's do not take spaces after the backslash
  73  * (like the c-shell), which will give you a very weird compiler
  74  * message.
  75  */
  76
  77 typedef struct t_inputrec_strings
  78 {
  79     char tcgrps[STRLEN], tau_t[STRLEN], ref_t[STRLEN],
  80          acc[STRLEN], accgrps[STRLEN], freeze[STRLEN], frdim[STRLEN],
  81          energy[STRLEN], user1[STRLEN], user2[STRLEN], vcm[STRLEN], x_compressed_groups[STRLEN],
  82          couple_moltype[STRLEN], orirefitgrp[STRLEN], egptable[STRLEN], egpexcl[STRLEN],
  83          wall_atomtype[STRLEN], wall_density[STRLEN], deform[STRLEN], QMMM[STRLEN],
  84          imd_grp[STRLEN];
  85     char   fep_lambda[efptNR][STRLEN];
  86     char   lambda_weights[STRLEN];
  87     char **pull_grp;
  88     char **rot_grp;
  89     char   anneal[STRLEN], anneal_npoints[STRLEN],
  90            anneal_time[STRLEN], anneal_temp[STRLEN];
  91     char   QMmethod[STRLEN], QMbasis[STRLEN], QMcharge[STRLEN], QMmult[STRLEN],
  92            bSH[STRLEN], CASorbitals[STRLEN], CASelectrons[STRLEN], SAon[STRLEN],
  93            SAoff[STRLEN], SAsteps[STRLEN], bTS[STRLEN], bOPT[STRLEN];
  94     char efield_x[STRLEN], efield_xt[STRLEN], efield_y[STRLEN],
  95          efield_yt[STRLEN], efield_z[STRLEN], efield_zt[STRLEN];
  96
  97 } gmx_inputrec_strings;
  98
  99 static gmx_inputrec_strings *is = NULL;
 100
 101 void init_inputrec_strings()
 102 {
 103     if (is)
 104     {
 105         gmx_incons("Attempted to call init_inputrec_strings before calling done_inputrec_strings. Only one inputrec (i.e. .mdp file) can be parsed at a time.");
 106     }
 107     snew(is, 1);
 108 }
 109
 110 void done_inputrec_strings()
 111 {
 112     sfree(is);
 113     is = NULL;
 114 }
 115
 116 static char swapgrp[STRLEN], splitgrp0[STRLEN], splitgrp1[STRLEN], solgrp[STRLEN];
 117
 118 enum {
 119     egrptpALL,         /* All particles have to be a member of a group.     */
 120     egrptpALL_GENREST, /* A rest group with name is generated for particles *
 121                         * that are not part of any group.                   */
 122     egrptpPART,        /* As egrptpALL_GENREST, but no name is generated    *
 123                         * for the rest group.                               */
 124     egrptpONE          /* Merge all selected groups into one group,         *
 125                         * make a rest group for the remaining particles.    */
 126 };
 127
 128 static const char *constraints[eshNR+1]    = {
 129     "none", "h-bonds", "all-bonds", "h-angles", "all-angles", NULL
 130 };
 131
 132 static const char *couple_lam[ecouplamNR+1]    = {
 133     "vdw-q", "vdw", "q", "none", NULL
 134 };
 135
 136 void init_ir(t_inputrec *ir, t_gromppopts *opts)
 137 {
 138     snew(opts->include, STRLEN);
 139     snew(opts->define, STRLEN);
 140     snew(ir->fepvals, 1);
 141     snew(ir->expandedvals, 1);
 142     snew(ir->simtempvals, 1);
 143 }
 144
 145 static void GetSimTemps(int ntemps, t_simtemp *simtemp, double *temperature_lambdas)
 146 {
 147
 148     int i;
 149
 150     for (i = 0; i < ntemps; i++)
 151     {
 152         /* simple linear scaling -- allows more control */
 153         if (simtemp->eSimTempScale == esimtempLINEAR)
 154         {
 155             simtemp->temperatures[i] = simtemp->simtemp_low + (simtemp->simtemp_high-simtemp->simtemp_low)*temperature_lambdas[i];
 156         }
 157         else if (simtemp->eSimTempScale == esimtempGEOMETRIC)  /* should give roughly equal acceptance for constant heat capacity . . . */
 158         {
 159             simtemp->temperatures[i] = simtemp->simtemp_low * pow(simtemp->simtemp_high/simtemp->simtemp_low, (1.0*i)/(ntemps-1));
 160         }
 161         else if (simtemp->eSimTempScale == esimtempEXPONENTIAL)
 162         {
 163             simtemp->temperatures[i] = simtemp->simtemp_low + (simtemp->simtemp_high-simtemp->simtemp_low)*((exp(temperature_lambdas[i])-1)/(exp(1.0)-1));
 164         }
 165         else
 166         {
 167             char errorstr[128];
 168             sprintf(errorstr, "eSimTempScale=%d not defined", simtemp->eSimTempScale);
 169             gmx_fatal(FARGS, errorstr);
 170         }
 171     }
 172 }
 173
 174
 175
 176 static void _low_check(gmx_bool b, char *s, warninp_t wi)
 177 {
 178     if (b)
 179     {
 180         warning_error(wi, s);
 181     }
 182 }
 183
 184 static void check_nst(const char *desc_nst, int nst,
 185                       const char *desc_p, int *p,
 186                       warninp_t wi)
 187 {
 188     char buf[STRLEN];
 189
 190     if (*p > 0 && *p % nst != 0)
 191     {
 192         /* Round up to the next multiple of nst */
 193         *p = ((*p)/nst + 1)*nst;
 194         sprintf(buf, "%s should be a multiple of %s, changing %s to %d\n",
 195                 desc_p, desc_nst, desc_p, *p);
 196         warning(wi, buf);
 197     }
 198 }
 199
 200 static gmx_bool ir_NVE(const t_inputrec *ir)
 201 {
 202     return ((ir->eI == eiMD || EI_VV(ir->eI)) && ir->etc == etcNO);
 203 }
 204
 205 static int lcd(int n1, int n2)
 206 {
 207     int d, i;
 208
 209     d = 1;
 210     for (i = 2; (i <= n1 && i <= n2); i++)
 211     {
 212         if (n1 % i == 0 && n2 % i == 0)
 213         {
 214             d = i;
 215         }
 216     }
 217
 218     return d;
 219 }
 220
 221 static void process_interaction_modifier(const t_inputrec *ir, int *eintmod)
 222 {
 223     if (*eintmod == eintmodPOTSHIFT_VERLET)
 224     {
 225         if (ir->cutoff_scheme == ecutsVERLET)
 226         {
 227             *eintmod = eintmodPOTSHIFT;
 228         }
 229         else
 230         {
 231             *eintmod = eintmodNONE;
 232         }
 233     }
 234 }
 235
 236 void check_ir(const char *mdparin, t_inputrec *ir, t_gromppopts *opts,
 237               warninp_t wi)
 238 /* Check internal consistency.
 239  * NOTE: index groups are not set here yet, don't check things
 240  * like temperature coupling group options here, but in triple_check
 241  */
 242 {
 243     /* Strange macro: first one fills the err_buf, and then one can check
 244      * the condition, which will print the message and increase the error
 245      * counter.
 246      */
 247 #define CHECK(b) _low_check(b, err_buf, wi)
 248     char        err_buf[256], warn_buf[STRLEN];
 249     int         i, j;
 250     int         ns_type  = 0;
 251     real        dt_coupl = 0;
 252     real        dt_pcoupl;
 253     int         nstcmin;
 254     t_lambda   *fep    = ir->fepvals;
 255     t_expanded *expand = ir->expandedvals;
 256
 257     set_warning_line(wi, mdparin, -1);
 258
 259     /* BASIC CUT-OFF STUFF */
 260     if (ir->rcoulomb < 0)
 261     {
 262         warning_error(wi, "rcoulomb should be >= 0");
 263     }
 264     if (ir->rvdw < 0)
 265     {
 266         warning_error(wi, "rvdw should be >= 0");
 267     }
 268     if (ir->rlist < 0 &&
 269         !(ir->cutoff_scheme == ecutsVERLET && ir->verletbuf_tol > 0))
 270     {
 271         warning_error(wi, "rlist should be >= 0");
 272     }
 273
 274     process_interaction_modifier(ir, &ir->coulomb_modifier);
 275     process_interaction_modifier(ir, &ir->vdw_modifier);
 276
 277     if (ir->cutoff_scheme == ecutsGROUP)
 278     {
 279         warning_note(wi,
 280                      "The group cutoff scheme is deprecated in Gromacs 5.0 and will be removed in a future "
 281                      "release when all interaction forms are supported for the verlet scheme. The verlet "
 282                      "scheme already scales better, and it is compatible with GPUs and other accelerators.");
 283
 284         /* BASIC CUT-OFF STUFF */
 285         if (ir->rlist == 0 ||
 286             !((ir_coulomb_might_be_zero_at_cutoff(ir) && ir->rcoulomb > ir->rlist) ||
 287               (ir_vdw_might_be_zero_at_cutoff(ir)     && ir->rvdw     > ir->rlist)))
 288         {
 289             /* No switched potential and/or no twin-range:
 290              * we can set the long-range cut-off to the maximum of the other cut-offs.
 291              */
 292             ir->rlistlong = max_cutoff(ir->rlist, max_cutoff(ir->rvdw, ir->rcoulomb));
 293         }
 294         else if (ir->rlistlong < 0)
 295         {
 296             ir->rlistlong = max_cutoff(ir->rlist, max_cutoff(ir->rvdw, ir->rcoulomb));
 297             sprintf(warn_buf, "rlistlong was not set, setting it to %g (no buffer)",
 298                     ir->rlistlong);
 299             warning(wi, warn_buf);
 300         }
 301         if (ir->rlistlong == 0 && ir->ePBC != epbcNONE)
 302         {
 303             warning_error(wi, "Can not have an infinite cut-off with PBC");
 304         }
 305         if (ir->rlistlong > 0 && (ir->rlist == 0 || ir->rlistlong < ir->rlist))
 306         {
 307             warning_error(wi, "rlistlong can not be shorter than rlist");
 308         }
 309         if (IR_TWINRANGE(*ir) && ir->nstlist <= 0)
 310         {
 311             warning_error(wi, "Can not have nstlist<=0 with twin-range interactions");
 312         }
 313     }
 314
 315     if (ir->rlistlong == ir->rlist)
 316     {
 317         ir->nstcalclr = 0;
 318     }
 319     else if (ir->rlistlong > ir->rlist && ir->nstcalclr == 0)
 320     {
 321         warning_error(wi, "With different cutoffs for electrostatics and VdW, nstcalclr must be -1 or a positive number");
 322     }
 323
 324     if (ir->cutoff_scheme == ecutsVERLET)
 325     {
 326         real rc_max;
 327
 328         /* Normal Verlet type neighbor-list, currently only limited feature support */
 329         if (inputrec2nboundeddim(ir) < 3)
 330         {
 331             warning_error(wi, "With Verlet lists only full pbc or pbc=xy with walls is supported");
 332         }
 333         if (ir->rcoulomb != ir->rvdw)
 334         {
 335             warning_error(wi, "With Verlet lists rcoulomb!=rvdw is not supported");
 336         }
 337         if (ir->vdwtype == evdwSHIFT || ir->vdwtype == evdwSWITCH)
 338         {
 339             if (ir->vdw_modifier == eintmodNONE ||
 340                 ir->vdw_modifier == eintmodPOTSHIFT)
 341             {
 342                 ir->vdw_modifier = (ir->vdwtype == evdwSHIFT ? eintmodFORCESWITCH : eintmodPOTSWITCH);
 343
 344                 sprintf(warn_buf, "Replacing vdwtype=%s by the equivalent combination of vdwtype=%s and vdw_modifier=%s", evdw_names[ir->vdwtype], evdw_names[evdwCUT], eintmod_names[ir->vdw_modifier]);
 345                 warning_note(wi, warn_buf);
 346
 347                 ir->vdwtype = evdwCUT;
 348             }
 349             else
 350             {
 351                 sprintf(warn_buf, "Unsupported combination of vdwtype=%s and vdw_modifier=%s", evdw_names[ir->vdwtype], eintmod_names[ir->vdw_modifier]);
 352                 warning_error(wi, warn_buf);
 353             }
 354         }
 355
 356         if (!(ir->vdwtype == evdwCUT || ir->vdwtype == evdwPME))
 357         {
 358             warning_error(wi, "With Verlet lists only cut-off and PME LJ interactions are supported");
 359         }
 360         if (!(ir->coulombtype == eelCUT ||
 361               (EEL_RF(ir->coulombtype) && ir->coulombtype != eelRF_NEC) ||
 362               EEL_PME(ir->coulombtype) || ir->coulombtype == eelEWALD))
 363         {
 364             warning_error(wi, "With Verlet lists only cut-off, reaction-field, PME and Ewald electrostatics are supported");
 365         }
 366         if (!(ir->coulomb_modifier == eintmodNONE ||
 367               ir->coulomb_modifier == eintmodPOTSHIFT))
 368         {
 369             sprintf(warn_buf, "coulomb_modifier=%s is not supported with the Verlet cut-off scheme", eintmod_names[ir->coulomb_modifier]);
 370             warning_error(wi, warn_buf);
 371         }
 372
 373         if (ir->nstlist <= 0)
 374         {
 375             warning_error(wi, "With Verlet lists nstlist should be larger than 0");
 376         }
 377
 378         if (ir->nstlist < 10)
 379         {
 380             warning_note(wi, "With Verlet lists the optimal nstlist is >= 10, with GPUs >= 20. Note that with the Verlet scheme, nstlist has no effect on the accuracy of your simulation.");
 381         }
 382
 383         rc_max = max(ir->rvdw, ir->rcoulomb);
 384
 385         if (ir->verletbuf_tol <= 0)
 386         {
 387             if (ir->verletbuf_tol == 0)
 388             {
 389                 warning_error(wi, "Can not have Verlet buffer tolerance of exactly 0");
 390             }
 391
 392             if (ir->rlist < rc_max)
 393             {
 394                 warning_error(wi, "With verlet lists rlist can not be smaller than rvdw or rcoulomb");
 395             }
 396
 397             if (ir->rlist == rc_max && ir->nstlist > 1)
 398             {
 399                 warning_note(wi, "rlist is equal to rvdw and/or rcoulomb: there is no explicit Verlet buffer. The cluster pair list does have a buffering effect, but choosing a larger rlist might be necessary for good energy conservation.");
 400             }
 401         }
 402         else
 403         {
 404             if (ir->rlist > rc_max)
 405             {
 406                 warning_note(wi, "You have set rlist larger than the interaction cut-off, but you also have verlet-buffer-tolerance > 0. Will set rlist using verlet-buffer-tolerance.");
 407             }
 408
 409             if (ir->nstlist == 1)
 410             {
 411                 /* No buffer required */
 412                 ir->rlist = rc_max;
 413             }
 414             else
 415             {
 416                 if (EI_DYNAMICS(ir->eI))
 417                 {
 418                     if (inputrec2nboundeddim(ir) < 3)
 419                     {
 420                         warning_error(wi, "The box volume is required for calculating rlist from the energy drift with verlet-buffer-tolerance > 0. You are using at least one unbounded dimension, so no volume can be computed. Either use a finite box, or set rlist yourself together with verlet-buffer-tolerance = -1.");
 421                     }
 422                     /* Set rlist temporarily so we can continue processing */
 423                     ir->rlist = rc_max;
 424                 }
 425                 else
 426                 {
 427                     /* Set the buffer to 5% of the cut-off */
 428                     ir->rlist = (1.0 + verlet_buffer_ratio_nodynamics)*rc_max;
 429                 }
 430             }
 431         }
 432
 433         /* No twin-range calculations with Verlet lists */
 434         ir->rlistlong = ir->rlist;
 435     }
 436
 437     if (ir->nstcalclr == -1)
 438     {
 439         /* if rlist=rlistlong, this will later be changed to nstcalclr=0 */
 440         ir->nstcalclr = ir->nstlist;
 441     }
 442     else if (ir->nstcalclr > 0)
 443     {
 444         if (ir->nstlist > 0 && (ir->nstlist % ir->nstcalclr != 0))
 445         {
 446             warning_error(wi, "nstlist must be evenly divisible by nstcalclr. Use nstcalclr = -1 to automatically follow nstlist");
 447         }
 448     }
 449     else if (ir->nstcalclr < -1)
 450     {
 451         warning_error(wi, "nstcalclr must be a positive number (divisor of nstcalclr), or -1 to follow nstlist.");
 452     }
 453
 454     if (EEL_PME(ir->coulombtype) && ir->rcoulomb > ir->rvdw && ir->nstcalclr > 1)
 455     {
 456         warning_error(wi, "When used with PME, the long-range component of twin-range interactions must be updated every step (nstcalclr)");
 457     }
 458
 459     /* GENERAL INTEGRATOR STUFF */
 460     if (!(ir->eI == eiMD || EI_VV(ir->eI)))
 461     {
 462         ir->etc = etcNO;
 463     }
 464     if (ir->eI == eiVVAK)
 465     {
 466         sprintf(warn_buf, "Integrator method %s is implemented primarily for validation purposes; for molecular dynamics, you should probably be using %s or %s", ei_names[eiVVAK], ei_names[eiMD], ei_names[eiVV]);
 467         warning_note(wi, warn_buf);
 468     }
 469     if (!EI_DYNAMICS(ir->eI))
 470     {
 471         ir->epc = epcNO;
 472     }
 473     if (EI_DYNAMICS(ir->eI))
 474     {
 475         if (ir->nstcalcenergy < 0)
 476         {
 477             ir->nstcalcenergy = ir_optimal_nstcalcenergy(ir);
 478             if (ir->nstenergy != 0 && ir->nstenergy < ir->nstcalcenergy)
 479             {
 480                 /* nstcalcenergy larger than nstener does not make sense.
 481                  * We ideally want nstcalcenergy=nstener.
 482                  */
 483                 if (ir->nstlist > 0)
 484                 {
 485                     ir->nstcalcenergy = lcd(ir->nstenergy, ir->nstlist);
 486                 }
 487                 else
 488                 {
 489                     ir->nstcalcenergy = ir->nstenergy;
 490                 }
 491             }
 492         }
 493         else if ( (ir->nstenergy > 0 && ir->nstcalcenergy > ir->nstenergy) ||
 494                   (ir->efep != efepNO && ir->fepvals->nstdhdl > 0 &&
 495                    (ir->nstcalcenergy > ir->fepvals->nstdhdl) ) )
 496
 497         {
 498             const char *nsten    = "nstenergy";
 499             const char *nstdh    = "nstdhdl";
 500             const char *min_name = nsten;
 501             int         min_nst  = ir->nstenergy;
 502
 503             /* find the smallest of ( nstenergy, nstdhdl ) */
 504             if (ir->efep != efepNO && ir->fepvals->nstdhdl > 0 &&
 505                 (ir->nstenergy == 0 || ir->fepvals->nstdhdl < ir->nstenergy))
 506             {
 507                 min_nst  = ir->fepvals->nstdhdl;
 508                 min_name = nstdh;
 509             }
 510             /* If the user sets nstenergy small, we should respect that */
 511             sprintf(warn_buf,
 512                     "Setting nstcalcenergy (%d) equal to %s (%d)",
 513                     ir->nstcalcenergy, min_name, min_nst);
 514             warning_note(wi, warn_buf);
 515             ir->nstcalcenergy = min_nst;
 516         }
 517
 518         if (ir->epc != epcNO)
 519         {
 520             if (ir->nstpcouple < 0)
 521             {
 522                 ir->nstpcouple = ir_optimal_nstpcouple(ir);
 523             }
 524         }
 525         if (IR_TWINRANGE(*ir))
 526         {
 527             check_nst("nstlist", ir->nstlist,
 528                       "nstcalcenergy", &ir->nstcalcenergy, wi);
 529             if (ir->epc != epcNO)
 530             {
 531                 check_nst("nstlist", ir->nstlist,
 532                           "nstpcouple", &ir->nstpcouple, wi);
 533             }
 534         }
 535
 536         if (ir->nstcalcenergy > 0)
 537         {
 538             if (ir->efep != efepNO)
 539             {
 540                 /* nstdhdl should be a multiple of nstcalcenergy */
 541                 check_nst("nstcalcenergy", ir->nstcalcenergy,
 542                           "nstdhdl", &ir->fepvals->nstdhdl, wi);
 543                 /* nstexpanded should be a multiple of nstcalcenergy */
 544                 check_nst("nstcalcenergy", ir->nstcalcenergy,
 545                           "nstexpanded", &ir->expandedvals->nstexpanded, wi);
 546             }
 547             /* for storing exact averages nstenergy should be
 548              * a multiple of nstcalcenergy
 549              */
 550             check_nst("nstcalcenergy", ir->nstcalcenergy,
 551                       "nstenergy", &ir->nstenergy, wi);
 552         }
 553     }
 554
 555     if (ir->nsteps == 0 && !ir->bContinuation)
 556     {
 557         warning_note(wi, "For a correct single-point energy evaluation with nsteps = 0, use continuation = yes to avoid constraining the input coordinates.");
 558     }
 559
 560     /* LD STUFF */
 561     if ((EI_SD(ir->eI) || ir->eI == eiBD) &&
 562         ir->bContinuation && ir->ld_seed != -1)
 563     {
 564         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)");
 565     }
 566
 567     /* TPI STUFF */
 568     if (EI_TPI(ir->eI))
 569     {
 570         sprintf(err_buf, "TPI only works with pbc = %s", epbc_names[epbcXYZ]);
 571         CHECK(ir->ePBC != epbcXYZ);
 572         sprintf(err_buf, "TPI only works with ns = %s", ens_names[ensGRID]);
 573         CHECK(ir->ns_type != ensGRID);
 574         sprintf(err_buf, "with TPI nstlist should be larger than zero");
 575         CHECK(ir->nstlist <= 0);
 576         sprintf(err_buf, "TPI does not work with full electrostatics other than PME");
 577         CHECK(EEL_FULL(ir->coulombtype) && !EEL_PME(ir->coulombtype));
 578     }
 579
 580     /* SHAKE / LINCS */
 581     if ( (opts->nshake > 0) && (opts->bMorse) )
 582     {
 583         sprintf(warn_buf,
 584                 "Using morse bond-potentials while constraining bonds is useless");
 585         warning(wi, warn_buf);
 586     }
 587
 588     if ((EI_SD(ir->eI) || ir->eI == eiBD) &&
 589         ir->bContinuation && ir->ld_seed != -1)
 590     {
 591         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)");
 592     }
 593     /* verify simulated tempering options */
 594
 595     if (ir->bSimTemp)
 596     {
 597         gmx_bool bAllTempZero = TRUE;
 598         for (i = 0; i < fep->n_lambda; i++)
 599         {
 600             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]);
 601             CHECK((fep->all_lambda[efptTEMPERATURE][i] < 0) || (fep->all_lambda[efptTEMPERATURE][i] > 1));
 602             if (fep->all_lambda[efptTEMPERATURE][i] > 0)
 603             {
 604                 bAllTempZero = FALSE;
 605             }
 606         }
 607         sprintf(err_buf, "if simulated tempering is on, temperature-lambdas may not be all zero");
 608         CHECK(bAllTempZero == TRUE);
 609
 610         sprintf(err_buf, "Simulated tempering is currently only compatible with md-vv");
 611         CHECK(ir->eI != eiVV);
 612
 613         /* check compatability of the temperature coupling with simulated tempering */
 614
 615         if (ir->etc == etcNOSEHOOVER)
 616         {
 617             sprintf(warn_buf, "Nose-Hoover based temperature control such as [%s] my not be entirelyconsistent with simulated tempering", etcoupl_names[ir->etc]);
 618             warning_note(wi, warn_buf);
 619         }
 620
 621         /* check that the temperatures make sense */
 622
 623         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);
 624         CHECK(ir->simtempvals->simtemp_high <= ir->simtempvals->simtemp_low);
 625
 626         sprintf(err_buf, "Higher simulated tempering temperature (%g) must be >= zero", ir->simtempvals->simtemp_high);
 627         CHECK(ir->simtempvals->simtemp_high <= 0);
 628
 629         sprintf(err_buf, "Lower simulated tempering temperature (%g) must be >= zero", ir->simtempvals->simtemp_low);
 630         CHECK(ir->simtempvals->simtemp_low <= 0);
 631     }
 632
 633     /* verify free energy options */
 634
 635     if (ir->efep != efepNO)
 636     {
 637         fep = ir->fepvals;
 638         sprintf(err_buf, "The soft-core power is %d and can only be 1 or 2",
 639                 fep->sc_power);
 640         CHECK(fep->sc_alpha != 0 && fep->sc_power != 1 && fep->sc_power != 2);
 641
 642         sprintf(err_buf, "The soft-core sc-r-power is %d and can only be 6 or 48",
 643                 (int)fep->sc_r_power);
 644         CHECK(fep->sc_alpha != 0 && fep->sc_r_power != 6.0 && fep->sc_r_power != 48.0);
 645
 646         sprintf(err_buf, "Can't use postive delta-lambda (%g) if initial state/lambda does not start at zero", fep->delta_lambda);
 647         CHECK(fep->delta_lambda > 0 && ((fep->init_fep_state > 0) ||  (fep->init_lambda > 0)));
 648
 649         sprintf(err_buf, "Can't use postive delta-lambda (%g) with expanded ensemble simulations", fep->delta_lambda);
 650         CHECK(fep->delta_lambda > 0 && (ir->efep == efepEXPANDED));
 651
 652         sprintf(err_buf, "Can only use expanded ensemble with md-vv for now; should be supported for other integrators in 5.0");
 653         CHECK(!(EI_VV(ir->eI)) && (ir->efep == efepEXPANDED));
 654
 655         sprintf(err_buf, "Free-energy not implemented for Ewald");
 656         CHECK(ir->coulombtype == eelEWALD);
 657
 658         /* check validty of lambda inputs */
 659         if (fep->n_lambda == 0)
 660         {
 661             /* Clear output in case of no states:*/
 662             sprintf(err_buf, "init-lambda-state set to %d: no lambda states are defined.", fep->init_fep_state);
 663             CHECK((fep->init_fep_state >= 0) && (fep->n_lambda == 0));
 664         }
 665         else
 666         {
 667             sprintf(err_buf, "initial thermodynamic state %d does not exist, only goes to %d", fep->init_fep_state, fep->n_lambda-1);
 668             CHECK((fep->init_fep_state >= fep->n_lambda));
 669         }
 670
 671         sprintf(err_buf, "Lambda state must be set, either with init-lambda-state or with init-lambda");
 672         CHECK((fep->init_fep_state < 0) && (fep->init_lambda < 0));
 673
 674         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",
 675                 fep->init_lambda, fep->init_fep_state);
 676         CHECK((fep->init_fep_state >= 0) && (fep->init_lambda >= 0));
 677
 678
 679
 680         if ((fep->init_lambda >= 0) && (fep->delta_lambda == 0))
 681         {
 682             int n_lambda_terms;
 683             n_lambda_terms = 0;
 684             for (i = 0; i < efptNR; i++)
 685             {
 686                 if (fep->separate_dvdl[i])
 687                 {
 688                     n_lambda_terms++;
 689                 }
 690             }
 691             if (n_lambda_terms > 1)
 692             {
 693                 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.");
 694                 warning(wi, warn_buf);
 695             }
 696
 697             if (n_lambda_terms < 2 && fep->n_lambda > 0)
 698             {
 699                 warning_note(wi,
 700                              "init-lambda is deprecated for setting lambda state (except for slow growth). Use init-lambda-state instead.");
 701             }
 702         }
 703
 704         for (j = 0; j < efptNR; j++)
 705         {
 706             for (i = 0; i < fep->n_lambda; i++)
 707             {
 708                 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]);
 709                 CHECK((fep->all_lambda[j][i] < 0) || (fep->all_lambda[j][i] > 1));
 710             }
 711         }
 712
 713         if ((fep->sc_alpha > 0) && (!fep->bScCoul))
 714         {
 715             for (i = 0; i < fep->n_lambda; i++)
 716             {
 717                 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],
 718                         fep->all_lambda[efptCOUL][i]);
 719                 CHECK((fep->sc_alpha > 0) &&
 720                       (((fep->all_lambda[efptCOUL][i] > 0.0) &&
 721                         (fep->all_lambda[efptCOUL][i] < 1.0)) &&
 722                        ((fep->all_lambda[efptVDW][i] > 0.0) &&
 723                         (fep->all_lambda[efptVDW][i] < 1.0))));
 724             }
 725         }
 726
 727         if ((fep->bScCoul) && (EEL_PME(ir->coulombtype)))
 728         {
 729             real sigma, lambda, r_sc;
 730
 731             sigma  = 0.34;
 732             /* Maximum estimate for A and B charges equal with lambda power 1 */
 733             lambda = 0.5;
 734             r_sc   = pow(lambda*fep->sc_alpha*pow(sigma/ir->rcoulomb, fep->sc_r_power) + 1.0, 1.0/fep->sc_r_power);
 735             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.",
 736                     fep->sc_r_power,
 737                     sigma, lambda, r_sc - 1.0, ir->ewald_rtol);
 738             warning_note(wi, warn_buf);
 739         }
 740
 741         /*  Free Energy Checks -- In an ideal world, slow growth and FEP would
 742             be treated differently, but that's the next step */
 743
 744         for (i = 0; i < efptNR; i++)
 745         {
 746             for (j = 0; j < fep->n_lambda; j++)
 747             {
 748                 sprintf(err_buf, "%s[%d] must be between 0 and 1", efpt_names[i], j);
 749                 CHECK((fep->all_lambda[i][j] < 0) || (fep->all_lambda[i][j] > 1));
 750             }
 751         }
 752     }
 753
 754     if ((ir->bSimTemp) || (ir->efep == efepEXPANDED))
 755     {
 756         fep    = ir->fepvals;
 757         expand = ir->expandedvals;
 758
 759         /* checking equilibration of weights inputs for validity */
 760
 761         sprintf(err_buf, "weight-equil-number-all-lambda (%d) is ignored if lmc-weights-equil is not equal to %s",
 762                 expand->equil_n_at_lam, elmceq_names[elmceqNUMATLAM]);
 763         CHECK((expand->equil_n_at_lam > 0) && (expand->elmceq != elmceqNUMATLAM));
 764
 765         sprintf(err_buf, "weight-equil-number-samples (%d) is ignored if lmc-weights-equil is not equal to %s",
 766                 expand->equil_samples, elmceq_names[elmceqSAMPLES]);
 767         CHECK((expand->equil_samples > 0) && (expand->elmceq != elmceqSAMPLES));
 768
 769         sprintf(err_buf, "weight-equil-number-steps (%d) is ignored if lmc-weights-equil is not equal to %s",
 770                 expand->equil_steps, elmceq_names[elmceqSTEPS]);
 771         CHECK((expand->equil_steps > 0) && (expand->elmceq != elmceqSTEPS));
 772
 773         sprintf(err_buf, "weight-equil-wl-delta (%d) is ignored if lmc-weights-equil is not equal to %s",
 774                 expand->equil_samples, elmceq_names[elmceqWLDELTA]);
 775         CHECK((expand->equil_wl_delta > 0) && (expand->elmceq != elmceqWLDELTA));
 776
 777         sprintf(err_buf, "weight-equil-count-ratio (%f) is ignored if lmc-weights-equil is not equal to %s",
 778                 expand->equil_ratio, elmceq_names[elmceqRATIO]);
 779         CHECK((expand->equil_ratio > 0) && (expand->elmceq != elmceqRATIO));
 780
 781         sprintf(err_buf, "weight-equil-number-all-lambda (%d) must be a positive integer if lmc-weights-equil=%s",
 782                 expand->equil_n_at_lam, elmceq_names[elmceqNUMATLAM]);
 783         CHECK((expand->equil_n_at_lam <= 0) && (expand->elmceq == elmceqNUMATLAM));
 784
 785         sprintf(err_buf, "weight-equil-number-samples (%d) must be a positive integer if lmc-weights-equil=%s",
 786                 expand->equil_samples, elmceq_names[elmceqSAMPLES]);
 787         CHECK((expand->equil_samples <= 0) && (expand->elmceq == elmceqSAMPLES));
 788
 789         sprintf(err_buf, "weight-equil-number-steps (%d) must be a positive integer if lmc-weights-equil=%s",
 790                 expand->equil_steps, elmceq_names[elmceqSTEPS]);
 791         CHECK((expand->equil_steps <= 0) && (expand->elmceq == elmceqSTEPS));
 792
 793         sprintf(err_buf, "weight-equil-wl-delta (%f) must be > 0 if lmc-weights-equil=%s",
 794                 expand->equil_wl_delta, elmceq_names[elmceqWLDELTA]);
 795         CHECK((expand->equil_wl_delta <= 0) && (expand->elmceq == elmceqWLDELTA));
 796
 797         sprintf(err_buf, "weight-equil-count-ratio (%f) must be > 0 if lmc-weights-equil=%s",
 798                 expand->equil_ratio, elmceq_names[elmceqRATIO]);
 799         CHECK((expand->equil_ratio <= 0) && (expand->elmceq == elmceqRATIO));
 800
 801         sprintf(err_buf, "lmc-weights-equil=%s only possible when lmc-stats = %s or lmc-stats %s",
 802                 elmceq_names[elmceqWLDELTA], elamstats_names[elamstatsWL], elamstats_names[elamstatsWWL]);
 803         CHECK((expand->elmceq == elmceqWLDELTA) && (!EWL(expand->elamstats)));
 804
 805         sprintf(err_buf, "lmc-repeats (%d) must be greater than 0", expand->lmc_repeats);
 806         CHECK((expand->lmc_repeats <= 0));
 807         sprintf(err_buf, "minimum-var-min (%d) must be greater than 0", expand->minvarmin);
 808         CHECK((expand->minvarmin <= 0));
 809         sprintf(err_buf, "weight-c-range (%d) must be greater or equal to 0", expand->c_range);
 810         CHECK((expand->c_range < 0));
 811         sprintf(err_buf, "init-lambda-state (%d) must be zero if lmc-forced-nstart (%d)> 0 and lmc-move != 'no'",
 812                 fep->init_fep_state, expand->lmc_forced_nstart);
 813         CHECK((fep->init_fep_state != 0) && (expand->lmc_forced_nstart > 0) && (expand->elmcmove != elmcmoveNO));
 814         sprintf(err_buf, "lmc-forced-nstart (%d) must not be negative", expand->lmc_forced_nstart);
 815         CHECK((expand->lmc_forced_nstart < 0));
 816         sprintf(err_buf, "init-lambda-state (%d) must be in the interval [0,number of lambdas)", fep->init_fep_state);
 817         CHECK((fep->init_fep_state < 0) || (fep->init_fep_state >= fep->n_lambda));
 818
 819         sprintf(err_buf, "init-wl-delta (%f) must be greater than or equal to 0", expand->init_wl_delta);
 820         CHECK((expand->init_wl_delta < 0));
 821         sprintf(err_buf, "wl-ratio (%f) must be between 0 and 1", expand->wl_ratio);
 822         CHECK((expand->wl_ratio <= 0) || (expand->wl_ratio >= 1));
 823         sprintf(err_buf, "wl-scale (%f) must be between 0 and 1", expand->wl_scale);
 824         CHECK((expand->wl_scale <= 0) || (expand->wl_scale >= 1));
 825
 826         /* if there is no temperature control, we need to specify an MC temperature */
 827         sprintf(err_buf, "If there is no temperature control, and lmc-mcmove!= 'no',mc_temperature must be set to a positive number");
 828         if (expand->nstTij > 0)
 829         {
 830             sprintf(err_buf, "nst-transition-matrix (%d) must be an integer multiple of nstlog (%d)",
 831                     expand->nstTij, ir->nstlog);
 832             CHECK((mod(expand->nstTij, ir->nstlog) != 0));
 833         }
 834     }
 835
 836     /* PBC/WALLS */
 837     sprintf(err_buf, "walls only work with pbc=%s", epbc_names[epbcXY]);
 838     CHECK(ir->nwall && ir->ePBC != epbcXY);
 839
 840     /* VACUUM STUFF */
 841     if (ir->ePBC != epbcXYZ && ir->nwall != 2)
 842     {
 843         if (ir->ePBC == epbcNONE)
 844         {
 845             if (ir->epc != epcNO)
 846             {
 847                 warning(wi, "Turning off pressure coupling for vacuum system");
 848                 ir->epc = epcNO;
 849             }
 850         }
 851         else
 852         {
 853             sprintf(err_buf, "Can not have pressure coupling with pbc=%s",
 854                     epbc_names[ir->ePBC]);
 855             CHECK(ir->epc != epcNO);
 856         }
 857         sprintf(err_buf, "Can not have Ewald with pbc=%s", epbc_names[ir->ePBC]);
 858         CHECK(EEL_FULL(ir->coulombtype));
 859
 860         sprintf(err_buf, "Can not have dispersion correction with pbc=%s",
 861                 epbc_names[ir->ePBC]);
 862         CHECK(ir->eDispCorr != edispcNO);
 863     }
 864
 865     if (ir->rlist == 0.0)
 866     {
 867         sprintf(err_buf, "can only have neighborlist cut-off zero (=infinite)\n"
 868                 "with coulombtype = %s or coulombtype = %s\n"
 869                 "without periodic boundary conditions (pbc = %s) and\n"
 870                 "rcoulomb and rvdw set to zero",
 871                 eel_names[eelCUT], eel_names[eelUSER], epbc_names[epbcNONE]);
 872         CHECK(((ir->coulombtype != eelCUT) && (ir->coulombtype != eelUSER)) ||
 873               (ir->ePBC     != epbcNONE) ||
 874               (ir->rcoulomb != 0.0)      || (ir->rvdw != 0.0));
 875
 876         if (ir->nstlist < 0)
 877         {
 878             warning_error(wi, "Can not have heuristic neighborlist updates without cut-off");
 879         }
 880         if (ir->nstlist > 0)
 881         {
 882             warning_note(wi, "Simulating without cut-offs can be (slightly) faster with nstlist=0, nstype=simple and only one MPI rank");
 883         }
 884     }
 885
 886     /* COMM STUFF */
 887     if (ir->nstcomm == 0)
 888     {
 889         ir->comm_mode = ecmNO;
 890     }
 891     if (ir->comm_mode != ecmNO)
 892     {
 893         if (ir->nstcomm < 0)
 894         {
 895             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");
 896             ir->nstcomm = abs(ir->nstcomm);
 897         }
 898
 899         if (ir->nstcalcenergy > 0 && ir->nstcomm < ir->nstcalcenergy)
 900         {
 901             warning_note(wi, "nstcomm < nstcalcenergy defeats the purpose of nstcalcenergy, setting nstcomm to nstcalcenergy");
 902             ir->nstcomm = ir->nstcalcenergy;
 903         }
 904
 905         if (ir->comm_mode == ecmANGULAR)
 906         {
 907             sprintf(err_buf, "Can not remove the rotation around the center of mass with periodic molecules");
 908             CHECK(ir->bPeriodicMols);
 909             if (ir->ePBC != epbcNONE)
 910             {
 911                 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).");
 912             }
 913         }
 914     }
 915
 916     if (EI_STATE_VELOCITY(ir->eI) && ir->ePBC == epbcNONE && ir->comm_mode != ecmANGULAR)
 917     {
 918         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.");
 919     }
 920
 921     sprintf(err_buf, "Twin-range neighbour searching (NS) with simple NS"
 922             " algorithm not implemented");
 923     CHECK(((ir->rcoulomb > ir->rlist) || (ir->rvdw > ir->rlist))
 924           && (ir->ns_type == ensSIMPLE));
 925
 926     /* TEMPERATURE COUPLING */
 927     if (ir->etc == etcYES)
 928     {
 929         ir->etc = etcBERENDSEN;
 930         warning_note(wi, "Old option for temperature coupling given: "
 931                      "changing \"yes\" to \"Berendsen\"\n");
 932     }
 933
 934     if ((ir->etc == etcNOSEHOOVER) || (ir->epc == epcMTTK))
 935     {
 936         if (ir->opts.nhchainlength < 1)
 937         {
 938             sprintf(warn_buf, "number of Nose-Hoover chains (currently %d) cannot be less than 1,reset to 1\n", ir->opts.nhchainlength);
 939             ir->opts.nhchainlength = 1;
 940             warning(wi, warn_buf);
 941         }
 942
 943         if (ir->etc == etcNOSEHOOVER && !EI_VV(ir->eI) && ir->opts.nhchainlength > 1)
 944         {
 945             warning_note(wi, "leapfrog does not yet support Nose-Hoover chains, nhchainlength reset to 1");
 946             ir->opts.nhchainlength = 1;
 947         }
 948     }
 949     else
 950     {
 951         ir->opts.nhchainlength = 0;
 952     }
 953
 954     if (ir->eI == eiVVAK)
 955     {
 956         sprintf(err_buf, "%s implemented primarily for validation, and requires nsttcouple = 1 and nstpcouple = 1.",
 957                 ei_names[eiVVAK]);
 958         CHECK((ir->nsttcouple != 1) || (ir->nstpcouple != 1));
 959     }
 960
 961     if (ETC_ANDERSEN(ir->etc))
 962     {
 963         sprintf(err_buf, "%s temperature control not supported for integrator %s.", etcoupl_names[ir->etc], ei_names[ir->eI]);
 964         CHECK(!(EI_VV(ir->eI)));
 965
 966         if (ir->nstcomm > 0 && (ir->etc == etcANDERSEN))
 967         {
 968             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]);
 969             warning_note(wi, warn_buf);
 970         }
 971
 972         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]);
 973         CHECK(ir->nstcomm > 1 && (ir->etc == etcANDERSEN));
 974     }
 975
 976     if (ir->etc == etcBERENDSEN)
 977     {
 978         sprintf(warn_buf, "The %s thermostat does not generate the correct kinetic energy distribution. You might want to consider using the %s thermostat.",
 979                 ETCOUPLTYPE(ir->etc), ETCOUPLTYPE(etcVRESCALE));
 980         warning_note(wi, warn_buf);
 981     }
 982
 983     if ((ir->etc == etcNOSEHOOVER || ETC_ANDERSEN(ir->etc))
 984         && ir->epc == epcBERENDSEN)
 985     {
 986         sprintf(warn_buf, "Using Berendsen pressure coupling invalidates the "
 987                 "true ensemble for the thermostat");
 988         warning(wi, warn_buf);
 989     }
 990
 991     /* PRESSURE COUPLING */
 992     if (ir->epc == epcISOTROPIC)
 993     {
 994         ir->epc = epcBERENDSEN;
 995         warning_note(wi, "Old option for pressure coupling given: "
 996                      "changing \"Isotropic\" to \"Berendsen\"\n");
 997     }
 998
 999     if (ir->epc != epcNO)
1000     {
1001         dt_pcoupl = ir->nstpcouple*ir->delta_t;
1002
1003         sprintf(err_buf, "tau-p must be > 0 instead of %g\n", ir->tau_p);
1004         CHECK(ir->tau_p <= 0);
1005
1006         if (ir->tau_p/dt_pcoupl < pcouple_min_integration_steps(ir->epc))
1007         {
1008             sprintf(warn_buf, "For proper integration of the %s barostat, tau-p (%g) should be at least %d times larger than nstpcouple*dt (%g)",
1009                     EPCOUPLTYPE(ir->epc), ir->tau_p, pcouple_min_integration_steps(ir->epc), dt_pcoupl);
1010             warning(wi, warn_buf);
1011         }
1012
1013         sprintf(err_buf, "compressibility must be > 0 when using pressure"
1014                 " coupling %s\n", EPCOUPLTYPE(ir->epc));
1015         CHECK(ir->compress[XX][XX] < 0 || ir->compress[YY][YY] < 0 ||
1016               ir->compress[ZZ][ZZ] < 0 ||
1017               (trace(ir->compress) == 0 && ir->compress[YY][XX] <= 0 &&
1018                ir->compress[ZZ][XX] <= 0 && ir->compress[ZZ][YY] <= 0));
1019
1020         if (epcPARRINELLORAHMAN == ir->epc && opts->bGenVel)
1021         {
1022             sprintf(warn_buf,
1023                     "You are generating velocities so I am assuming you "
1024                     "are equilibrating a system. You are using "
1025                     "%s pressure coupling, but this can be "
1026                     "unstable for equilibration. If your system crashes, try "
1027                     "equilibrating first with Berendsen pressure coupling. If "
1028                     "you are not equilibrating the system, you can probably "
1029                     "ignore this warning.",
1030                     epcoupl_names[ir->epc]);
1031             warning(wi, warn_buf);
1032         }
1033     }
1034
1035     if (EI_VV(ir->eI))
1036     {
1037         if (ir->epc > epcNO)
1038         {
1039             if ((ir->epc != epcBERENDSEN) && (ir->epc != epcMTTK))
1040             {
1041                 warning_error(wi, "for md-vv and md-vv-avek, can only use Berendsen and Martyna-Tuckerman-Tobias-Klein (MTTK) equations for pressure control; MTTK is equivalent to Parrinello-Rahman.");
1042             }
1043         }
1044     }
1045     else
1046     {
1047         if (ir->epc == epcMTTK)
1048         {
1049             warning_error(wi, "MTTK pressure coupling requires a Velocity-verlet integrator");
1050         }
1051     }
1052
1053     /* ELECTROSTATICS */
1054     /* More checks are in triple check (grompp.c) */
1055
1056     if (ir->coulombtype == eelSWITCH)
1057     {
1058         sprintf(warn_buf, "coulombtype = %s is only for testing purposes and can lead to serious "
1059                 "artifacts, advice: use coulombtype = %s",
1060                 eel_names[ir->coulombtype],
1061                 eel_names[eelRF_ZERO]);
1062         warning(wi, warn_buf);
1063     }
1064
1065     if (ir->epsilon_r != 1 && ir->implicit_solvent == eisGBSA)
1066     {
1067         sprintf(warn_buf, "epsilon-r = %g with GB implicit solvent, will use this value for inner dielectric", ir->epsilon_r);
1068         warning_note(wi, warn_buf);
1069     }
1070
1071     if (EEL_RF(ir->coulombtype) && ir->epsilon_rf == 1 && ir->epsilon_r != 1)
1072     {
1073         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);
1074         warning(wi, warn_buf);
1075         ir->epsilon_rf = ir->epsilon_r;
1076         ir->epsilon_r  = 1.0;
1077     }
1078
1079     if (getenv("GMX_DO_GALACTIC_DYNAMICS") == NULL)
1080     {
1081         sprintf(err_buf, "epsilon-r must be >= 0 instead of %g\n", ir->epsilon_r);
1082         CHECK(ir->epsilon_r < 0);
1083     }
1084
1085     if (EEL_RF(ir->coulombtype))
1086     {
1087         /* reaction field (at the cut-off) */
1088
1089         if (ir->coulombtype == eelRF_ZERO)
1090         {
1091             sprintf(warn_buf, "With coulombtype = %s, epsilon-rf must be 0, assuming you meant epsilon_rf=0",
1092                     eel_names[ir->coulombtype]);
1093             CHECK(ir->epsilon_rf != 0);
1094             ir->epsilon_rf = 0.0;
1095         }
1096
1097         sprintf(err_buf, "epsilon-rf must be >= epsilon-r");
1098         CHECK((ir->epsilon_rf < ir->epsilon_r && ir->epsilon_rf != 0) ||
1099               (ir->epsilon_r == 0));
1100         if (ir->epsilon_rf == ir->epsilon_r)
1101         {
1102             sprintf(warn_buf, "Using epsilon-rf = epsilon-r with %s does not make sense",
1103                     eel_names[ir->coulombtype]);
1104             warning(wi, warn_buf);
1105         }
1106     }
1107     /* Allow rlist>rcoulomb for tabulated long range stuff. This just
1108      * means the interaction is zero outside rcoulomb, but it helps to
1109      * provide accurate energy conservation.
1110      */
1111     if (ir_coulomb_might_be_zero_at_cutoff(ir))
1112     {
1113         if (ir_coulomb_switched(ir))
1114         {
1115             sprintf(err_buf,
1116                     "With coulombtype = %s rcoulomb_switch must be < rcoulomb. Or, better: Use the potential modifier options!",
1117                     eel_names[ir->coulombtype]);
1118             CHECK(ir->rcoulomb_switch >= ir->rcoulomb);
1119         }
1120     }
1121     else if (ir->coulombtype == eelCUT || EEL_RF(ir->coulombtype))
1122     {
1123         if (ir->cutoff_scheme == ecutsGROUP && ir->coulomb_modifier == eintmodNONE)
1124         {
1125             sprintf(err_buf, "With coulombtype = %s, rcoulomb should be >= rlist unless you use a potential modifier",
1126                     eel_names[ir->coulombtype]);
1127             CHECK(ir->rlist > ir->rcoulomb);
1128         }
1129     }
1130
1131     if (ir->coulombtype == eelSWITCH || ir->coulombtype == eelSHIFT ||
1132         ir->vdwtype == evdwSWITCH || ir->vdwtype == evdwSHIFT)
1133     {
1134         sprintf(warn_buf,
1135                 "The switch/shift interaction settings are just for compatibility; you will get better "
1136                 "performance from applying potential modifiers to your interactions!\n");
1137         warning_note(wi, warn_buf);
1138     }
1139
1140     if (ir->coulombtype == eelPMESWITCH || ir->coulomb_modifier == eintmodPOTSWITCH)
1141     {
1142         if (ir->rcoulomb_switch/ir->rcoulomb < 0.9499)
1143         {
1144             real percentage  = 100*(ir->rcoulomb-ir->rcoulomb_switch)/ir->rcoulomb;
1145             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.",
1146                     percentage, ir->rcoulomb_switch, ir->rcoulomb, ir->ewald_rtol);
1147             warning(wi, warn_buf);
1148         }
1149     }
1150
1151     if (ir->vdwtype == evdwSWITCH || ir->vdw_modifier == eintmodPOTSWITCH)
1152     {
1153         if (ir->rvdw_switch == 0)
1154         {
1155             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.");
1156             warning(wi, warn_buf);
1157         }
1158     }
1159
1160     if (EEL_FULL(ir->coulombtype))
1161     {
1162         if (ir->coulombtype == eelPMESWITCH || ir->coulombtype == eelPMEUSER ||
1163             ir->coulombtype == eelPMEUSERSWITCH)
1164         {
1165             sprintf(err_buf, "With coulombtype = %s, rcoulomb must be <= rlist",
1166                     eel_names[ir->coulombtype]);
1167             CHECK(ir->rcoulomb > ir->rlist);
1168         }
1169         else if (ir->cutoff_scheme == ecutsGROUP && ir->coulomb_modifier == eintmodNONE)
1170         {
1171             if (ir->coulombtype == eelPME || ir->coulombtype == eelP3M_AD)
1172             {
1173                 sprintf(err_buf,
1174                         "With coulombtype = %s (without modifier), rcoulomb must be equal to rlist,\n"
1175                         "or rlistlong if nstcalclr=1. For optimal energy conservation,consider using\n"
1176                         "a potential modifier.", eel_names[ir->coulombtype]);
1177                 if (ir->nstcalclr == 1)
1178                 {
1179                     CHECK(ir->rcoulomb != ir->rlist && ir->rcoulomb != ir->rlistlong);
1180                 }
1181                 else
1182                 {
1183                     CHECK(ir->rcoulomb != ir->rlist);
1184                 }
1185             }
1186         }
1187     }
1188
1189     if (EEL_PME(ir->coulombtype) || EVDW_PME(ir->vdwtype))
1190     {
1191         if (ir->pme_order < 3)
1192         {
1193             warning_error(wi, "pme-order can not be smaller than 3");
1194         }
1195     }
1196
1197     if (ir->nwall == 2 && EEL_FULL(ir->coulombtype))
1198     {
1199         if (ir->ewald_geometry == eewg3D)
1200         {
1201             sprintf(warn_buf, "With pbc=%s you should use ewald-geometry=%s",
1202                     epbc_names[ir->ePBC], eewg_names[eewg3DC]);
1203             warning(wi, warn_buf);
1204         }
1205         /* This check avoids extra pbc coding for exclusion corrections */
1206         sprintf(err_buf, "wall-ewald-zfac should be >= 2");
1207         CHECK(ir->wall_ewald_zfac < 2);
1208     }
1209
1210     if (ir_vdw_switched(ir))
1211     {
1212         sprintf(err_buf, "With switched vdw forces or potentials, rvdw-switch must be < rvdw");
1213         CHECK(ir->rvdw_switch >= ir->rvdw);
1214
1215         if (ir->rvdw_switch < 0.5*ir->rvdw)
1216         {
1217             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.",
1218                     ir->rvdw_switch, ir->rvdw);
1219             warning_note(wi, warn_buf);
1220         }
1221     }
1222     else if (ir->vdwtype == evdwCUT || ir->vdwtype == evdwPME)
1223     {
1224         if (ir->cutoff_scheme == ecutsGROUP && ir->vdw_modifier == eintmodNONE)
1225         {
1226             sprintf(err_buf, "With vdwtype = %s, rvdw must be >= rlist unless you use a potential modifier", evdw_names[ir->vdwtype]);
1227             CHECK(ir->rlist > ir->rvdw);
1228         }
1229     }
1230
1231     if (ir->vdwtype == evdwPME)
1232     {
1233         if (!(ir->vdw_modifier == eintmodNONE || ir->vdw_modifier == eintmodPOTSHIFT))
1234         {
1235             sprintf(err_buf, "With vdwtype = %s, the only supported modifiers are %s a\
1236 nd %s",
1237                     evdw_names[ir->vdwtype],
1238                     eintmod_names[eintmodPOTSHIFT],
1239                     eintmod_names[eintmodNONE]);
1240         }
1241     }
1242
1243     if (ir->cutoff_scheme == ecutsGROUP)
1244     {
1245         if (((ir->coulomb_modifier != eintmodNONE && ir->rcoulomb == ir->rlist) ||
1246              (ir->vdw_modifier != eintmodNONE && ir->rvdw == ir->rlist)) &&
1247             ir->nstlist != 1)
1248         {
1249             warning_note(wi, "With exact cut-offs, rlist should be "
1250                          "larger than rcoulomb and rvdw, so that there "
1251                          "is a buffer region for particle motion "
1252                          "between neighborsearch steps");
1253         }
1254
1255         if (ir_coulomb_is_zero_at_cutoff(ir) && ir->rlistlong <= ir->rcoulomb)
1256         {
1257             sprintf(warn_buf, "For energy conservation with switch/shift potentials, %s should be 0.1 to 0.3 nm larger than rcoulomb.",
1258                     IR_TWINRANGE(*ir) ? "rlistlong" : "rlist");
1259             warning_note(wi, warn_buf);
1260         }
1261         if (ir_vdw_switched(ir) && (ir->rlistlong <= ir->rvdw))
1262         {
1263             sprintf(warn_buf, "For energy conservation with switch/shift potentials, %s should be 0.1 to 0.3 nm larger than rvdw.",
1264                     IR_TWINRANGE(*ir) ? "rlistlong" : "rlist");
1265             warning_note(wi, warn_buf);
1266         }
1267     }
1268
1269     if (ir->vdwtype == evdwUSER && ir->eDispCorr != edispcNO)
1270     {
1271         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.");
1272     }
1273
1274     if (ir->nstlist == -1)
1275     {
1276         sprintf(err_buf, "With nstlist=-1 rvdw and rcoulomb should be smaller than rlist to account for diffusion and possibly charge-group radii");
1277         CHECK(ir->rvdw >= ir->rlist || ir->rcoulomb >= ir->rlist);
1278     }
1279     sprintf(err_buf, "nstlist can not be smaller than -1");
1280     CHECK(ir->nstlist < -1);
1281
1282     if (ir->eI == eiLBFGS && (ir->coulombtype == eelCUT || ir->vdwtype == evdwCUT)
1283         && ir->rvdw != 0)
1284     {
1285         warning(wi, "For efficient BFGS minimization, use switch/shift/pme instead of cut-off.");
1286     }
1287
1288     if (ir->eI == eiLBFGS && ir->nbfgscorr <= 0)
1289     {
1290         warning(wi, "Using L-BFGS with nbfgscorr<=0 just gets you steepest descent.");
1291     }
1292
1293     /* ENERGY CONSERVATION */
1294     if (ir_NVE(ir) && ir->cutoff_scheme == ecutsGROUP)
1295     {
1296         if (!ir_vdw_might_be_zero_at_cutoff(ir) && ir->rvdw > 0 && ir->vdw_modifier == eintmodNONE)
1297         {
1298             sprintf(warn_buf, "You are using a cut-off for VdW interactions with NVE, for good energy conservation use vdwtype = %s (possibly with DispCorr)",
1299                     evdw_names[evdwSHIFT]);
1300             warning_note(wi, warn_buf);
1301         }
1302         if (!ir_coulomb_might_be_zero_at_cutoff(ir) && ir->rcoulomb > 0)
1303         {
1304             sprintf(warn_buf, "You are using a cut-off for electrostatics with NVE, for good energy conservation use coulombtype = %s or %s",
1305                     eel_names[eelPMESWITCH], eel_names[eelRF_ZERO]);
1306             warning_note(wi, warn_buf);
1307         }
1308     }
1309
1310     if (EI_VV(ir->eI) && IR_TWINRANGE(*ir) && ir->nstlist > 1)
1311     {
1312         sprintf(warn_buf, "Twin-range multiple time stepping does not work with integrator %s.", ei_names[ir->eI]);
1313         warning_error(wi, warn_buf);
1314     }
1315
1316     /* IMPLICIT SOLVENT */
1317     if (ir->coulombtype == eelGB_NOTUSED)
1318     {
1319         ir->coulombtype      = eelCUT;
1320         ir->implicit_solvent = eisGBSA;
1321         fprintf(stderr, "Note: Old option for generalized born electrostatics given:\n"
1322                 "Changing coulombtype from \"generalized-born\" to \"cut-off\" and instead\n"
1323                 "setting implicit-solvent value to \"GBSA\" in input section.\n");
1324     }
1325
1326     if (ir->sa_algorithm == esaSTILL)
1327     {
1328         sprintf(err_buf, "Still SA algorithm not available yet, use %s or %s instead\n", esa_names[esaAPPROX], esa_names[esaNO]);
1329         CHECK(ir->sa_algorithm == esaSTILL);
1330     }
1331
1332     if (ir->implicit_solvent == eisGBSA)
1333     {
1334         sprintf(err_buf, "With GBSA implicit solvent, rgbradii must be equal to rlist.");
1335         CHECK(ir->rgbradii != ir->rlist);
1336
1337         if (ir->coulombtype != eelCUT)
1338         {
1339             sprintf(err_buf, "With GBSA, coulombtype must be equal to %s\n", eel_names[eelCUT]);
1340             CHECK(ir->coulombtype != eelCUT);
1341         }
1342         if (ir->vdwtype != evdwCUT)
1343         {
1344             sprintf(err_buf, "With GBSA, vdw-type must be equal to %s\n", evdw_names[evdwCUT]);
1345             CHECK(ir->vdwtype != evdwCUT);
1346         }
1347         if (ir->nstgbradii < 1)
1348         {
1349             sprintf(warn_buf, "Using GBSA with nstgbradii<1, setting nstgbradii=1");
1350             warning_note(wi, warn_buf);
1351             ir->nstgbradii = 1;
1352         }
1353         if (ir->sa_algorithm == esaNO)
1354         {
1355             sprintf(warn_buf, "No SA (non-polar) calculation requested together with GB. Are you sure this is what you want?\n");
1356             warning_note(wi, warn_buf);
1357         }
1358         if (ir->sa_surface_tension < 0 && ir->sa_algorithm != esaNO)
1359         {
1360             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");
1361             warning_note(wi, warn_buf);
1362
1363             if (ir->gb_algorithm == egbSTILL)
1364             {
1365                 ir->sa_surface_tension = 0.0049 * CAL2JOULE * 100;
1366             }
1367             else
1368             {
1369                 ir->sa_surface_tension = 0.0054 * CAL2JOULE * 100;
1370             }
1371         }
1372         if (ir->sa_surface_tension == 0 && ir->sa_algorithm != esaNO)
1373         {
1374             sprintf(err_buf, "Surface tension set to 0 while SA-calculation requested\n");
1375             CHECK(ir->sa_surface_tension == 0 && ir->sa_algorithm != esaNO);
1376         }
1377
1378     }
1379
1380     if (ir->bAdress)
1381     {
1382         if (ir->cutoff_scheme != ecutsGROUP)
1383         {
1384             warning_error(wi, "AdresS simulation supports only cutoff-scheme=group");
1385         }
1386         if (!EI_SD(ir->eI))
1387         {
1388             warning_error(wi, "AdresS simulation supports only stochastic dynamics");
1389         }
1390         if (ir->epc != epcNO)
1391         {
1392             warning_error(wi, "AdresS simulation does not support pressure coupling");
1393         }
1394         if (EEL_FULL(ir->coulombtype))
1395         {
1396             warning_error(wi, "AdresS simulation does not support long-range electrostatics");
1397         }
1398     }
1399 }
1400
1401 /* count the number of text elemets separated by whitespace in a string.
1402     str = the input string
1403     maxptr = the maximum number of allowed elements
1404     ptr = the output array of pointers to the first character of each element
1405     returns: the number of elements. */
1406 int str_nelem(const char *str, int maxptr, char *ptr[])
1407 {
1408     int   np = 0;
1409     char *copy0, *copy;
1410
1411     copy0 = strdup(str);
1412     copy  = copy0;
1413     ltrim(copy);
1414     while (*copy != '\0')
1415     {
1416         if (np >= maxptr)
1417         {
1418             gmx_fatal(FARGS, "Too many groups on line: '%s' (max is %d)",
1419                       str, maxptr);
1420         }
1421         if (ptr)
1422         {
1423             ptr[np] = copy;
1424         }
1425         np++;
1426         while ((*copy != '\0') && !isspace(*copy))
1427         {
1428             copy++;
1429         }
1430         if (*copy != '\0')
1431         {
1432             *copy = '\0';
1433             copy++;
1434         }
1435         ltrim(copy);
1436     }
1437     if (ptr == NULL)
1438     {
1439         sfree(copy0);
1440     }
1441
1442     return np;
1443 }
1444
1445 /* interpret a number of doubles from a string and put them in an array,
1446    after allocating space for them.
1447    str = the input string
1448    n = the (pre-allocated) number of doubles read
1449    r = the output array of doubles. */
1450 static void parse_n_real(char *str, int *n, real **r)
1451 {
1452     char *ptr[MAXPTR];
1453     int   i;
1454
1455     *n = str_nelem(str, MAXPTR, ptr);
1456
1457     snew(*r, *n);
1458     for (i = 0; i < *n; i++)
1459     {
1460         (*r)[i] = strtod(ptr[i], NULL);
1461     }
1462 }
1463
1464 static void do_fep_params(t_inputrec *ir, char fep_lambda[][STRLEN], char weights[STRLEN])
1465 {
1466
1467     int         i, j, max_n_lambda, nweights, nfep[efptNR];
1468     t_lambda   *fep    = ir->fepvals;
1469     t_expanded *expand = ir->expandedvals;
1470     real      **count_fep_lambdas;
1471     gmx_bool    bOneLambda = TRUE;
1472
1473     snew(count_fep_lambdas, efptNR);
1474
1475     /* FEP input processing */
1476     /* first, identify the number of lambda values for each type.
1477        All that are nonzero must have the same number */
1478
1479     for (i = 0; i < efptNR; i++)
1480     {
1481         parse_n_real(fep_lambda[i], &(nfep[i]), &(count_fep_lambdas[i]));
1482     }
1483
1484     /* now, determine the number of components.  All must be either zero, or equal. */
1485
1486     max_n_lambda = 0;
1487     for (i = 0; i < efptNR; i++)
1488     {
1489         if (nfep[i] > max_n_lambda)
1490         {
1491             max_n_lambda = nfep[i];  /* here's a nonzero one.  All of them
1492                                         must have the same number if its not zero.*/
1493             break;
1494         }
1495     }
1496
1497     for (i = 0; i < efptNR; i++)
1498     {
1499         if (nfep[i] == 0)
1500         {
1501             ir->fepvals->separate_dvdl[i] = FALSE;
1502         }
1503         else if (nfep[i] == max_n_lambda)
1504         {
1505             if (i != efptTEMPERATURE)  /* we treat this differently -- not really a reason to compute the derivative with
1506                                           respect to the temperature currently */
1507             {
1508                 ir->fepvals->separate_dvdl[i] = TRUE;
1509             }
1510         }
1511         else
1512         {
1513             gmx_fatal(FARGS, "Number of lambdas (%d) for FEP type %s not equal to number of other types (%d)",
1514                       nfep[i], efpt_names[i], max_n_lambda);
1515         }
1516     }
1517     /* we don't print out dhdl if the temperature is changing, since we can't correctly define dhdl in this case */
1518     ir->fepvals->separate_dvdl[efptTEMPERATURE] = FALSE;
1519
1520     /* the number of lambdas is the number we've read in, which is either zero
1521        or the same for all */
1522     fep->n_lambda = max_n_lambda;
1523
1524     /* allocate space for the array of lambda values */
1525     snew(fep->all_lambda, efptNR);
1526     /* if init_lambda is defined, we need to set lambda */
1527     if ((fep->init_lambda > 0) && (fep->n_lambda == 0))
1528     {
1529         ir->fepvals->separate_dvdl[efptFEP] = TRUE;
1530     }
1531     /* otherwise allocate the space for all of the lambdas, and transfer the data */
1532     for (i = 0; i < efptNR; i++)
1533     {
1534         snew(fep->all_lambda[i], fep->n_lambda);
1535         if (nfep[i] > 0)  /* if it's zero, then the count_fep_lambda arrays
1536                              are zero */
1537         {
1538             for (j = 0; j < fep->n_lambda; j++)
1539             {
1540                 fep->all_lambda[i][j] = (double)count_fep_lambdas[i][j];
1541             }
1542             sfree(count_fep_lambdas[i]);
1543         }
1544     }
1545     sfree(count_fep_lambdas);
1546
1547     /* "fep-vals" is either zero or the full number. If zero, we'll need to define fep-lambdas for internal
1548        bookkeeping -- for now, init_lambda */
1549
1550     if ((nfep[efptFEP] == 0) && (fep->init_lambda >= 0))
1551     {
1552         for (i = 0; i < fep->n_lambda; i++)
1553         {
1554             fep->all_lambda[efptFEP][i] = fep->init_lambda;
1555         }
1556     }
1557
1558     /* check to see if only a single component lambda is defined, and soft core is defined.
1559        In this case, turn on coulomb soft core */
1560
1561     if (max_n_lambda == 0)
1562     {
1563         bOneLambda = TRUE;
1564     }
1565     else
1566     {
1567         for (i = 0; i < efptNR; i++)
1568         {
1569             if ((nfep[i] != 0) && (i != efptFEP))
1570             {
1571                 bOneLambda = FALSE;
1572             }
1573         }
1574     }
1575     if ((bOneLambda) && (fep->sc_alpha > 0))
1576     {
1577         fep->bScCoul = TRUE;
1578     }
1579
1580     /* Fill in the others with the efptFEP if they are not explicitly
1581        specified (i.e. nfep[i] == 0).  This means if fep is not defined,
1582        they are all zero. */
1583
1584     for (i = 0; i < efptNR; i++)
1585     {
1586         if ((nfep[i] == 0) && (i != efptFEP))
1587         {
1588             for (j = 0; j < fep->n_lambda; j++)
1589             {
1590                 fep->all_lambda[i][j] = fep->all_lambda[efptFEP][j];
1591             }
1592         }
1593     }
1594
1595
1596     /* make it easier if sc_r_power = 48 by increasing it to the 4th power, to be in the right scale. */
1597     if (fep->sc_r_power == 48)
1598     {
1599         if (fep->sc_alpha > 0.1)
1600         {
1601             gmx_fatal(FARGS, "sc_alpha (%f) for sc_r_power = 48 should usually be between 0.001 and 0.004", fep->sc_alpha);
1602         }
1603     }
1604
1605     expand = ir->expandedvals;
1606     /* now read in the weights */
1607     parse_n_real(weights, &nweights, &(expand->init_lambda_weights));
1608     if (nweights == 0)
1609     {
1610         snew(expand->init_lambda_weights, fep->n_lambda); /* initialize to zero */
1611     }
1612     else if (nweights != fep->n_lambda)
1613     {
1614         gmx_fatal(FARGS, "Number of weights (%d) is not equal to number of lambda values (%d)",
1615                   nweights, fep->n_lambda);
1616     }
1617     if ((expand->nstexpanded < 0) && (ir->efep != efepNO))
1618     {
1619         expand->nstexpanded = fep->nstdhdl;
1620         /* if you don't specify nstexpanded when doing expanded ensemble free energy calcs, it is set to nstdhdl */
1621     }
1622     if ((expand->nstexpanded < 0) && ir->bSimTemp)
1623     {
1624         expand->nstexpanded = 2*(int)(ir->opts.tau_t[0]/ir->delta_t);
1625         /* if you don't specify nstexpanded when doing expanded ensemble simulated tempering, it is set to
1626            2*tau_t just to be careful so it's not to frequent  */
1627     }
1628 }
1629
1630
1631 static void do_simtemp_params(t_inputrec *ir)
1632 {
1633
1634     snew(ir->simtempvals->temperatures, ir->fepvals->n_lambda);
1635     GetSimTemps(ir->fepvals->n_lambda, ir->simtempvals, ir->fepvals->all_lambda[efptTEMPERATURE]);
1636
1637     return;
1638 }
1639
1640 static void do_wall_params(t_inputrec *ir,
1641                            char *wall_atomtype, char *wall_density,
1642                            t_gromppopts *opts)
1643 {
1644     int    nstr, i;
1645     char  *names[MAXPTR];
1646     double dbl;
1647
1648     opts->wall_atomtype[0] = NULL;
1649     opts->wall_atomtype[1] = NULL;
1650
1651     ir->wall_atomtype[0] = -1;
1652     ir->wall_atomtype[1] = -1;
1653     ir->wall_density[0]  = 0;
1654     ir->wall_density[1]  = 0;
1655
1656     if (ir->nwall > 0)
1657     {
1658         nstr = str_nelem(wall_atomtype, MAXPTR, names);
1659         if (nstr != ir->nwall)
1660         {
1661             gmx_fatal(FARGS, "Expected %d elements for wall_atomtype, found %d",
1662                       ir->nwall, nstr);
1663         }
1664         for (i = 0; i < ir->nwall; i++)
1665         {
1666             opts->wall_atomtype[i] = strdup(names[i]);
1667         }
1668
1669         if (ir->wall_type == ewt93 || ir->wall_type == ewt104)
1670         {
1671             nstr = str_nelem(wall_density, MAXPTR, names);
1672             if (nstr != ir->nwall)
1673             {
1674                 gmx_fatal(FARGS, "Expected %d elements for wall-density, found %d", ir->nwall, nstr);
1675             }
1676             for (i = 0; i < ir->nwall; i++)
1677             {
1678                 sscanf(names[i], "%lf", &dbl);
1679                 if (dbl <= 0)
1680                 {
1681                     gmx_fatal(FARGS, "wall-density[%d] = %f\n", i, dbl);
1682                 }
1683                 ir->wall_density[i] = dbl;
1684             }
1685         }
1686     }
1687 }
1688
1689 static void add_wall_energrps(gmx_groups_t *groups, int nwall, t_symtab *symtab)
1690 {
1691     int     i;
1692     t_grps *grps;
1693     char    str[STRLEN];
1694
1695     if (nwall > 0)
1696     {
1697         srenew(groups->grpname, groups->ngrpname+nwall);
1698         grps = &(groups->grps[egcENER]);
1699         srenew(grps->nm_ind, grps->nr+nwall);
1700         for (i = 0; i < nwall; i++)
1701         {
1702             sprintf(str, "wall%d", i);
1703             groups->grpname[groups->ngrpname] = put_symtab(symtab, str);
1704             grps->nm_ind[grps->nr++]          = groups->ngrpname++;
1705         }
1706     }
1707 }
1708
1709 void read_expandedparams(int *ninp_p, t_inpfile **inp_p,
1710                          t_expanded *expand, warninp_t wi)
1711 {
1712     int        ninp, nerror = 0;
1713     t_inpfile *inp;
1714
1715     ninp   = *ninp_p;
1716     inp    = *inp_p;
1717
1718     /* read expanded ensemble parameters */
1719     CCTYPE ("expanded ensemble variables");
1720     ITYPE ("nstexpanded", expand->nstexpanded, -1);
1721     EETYPE("lmc-stats", expand->elamstats, elamstats_names);
1722     EETYPE("lmc-move", expand->elmcmove, elmcmove_names);
1723     EETYPE("lmc-weights-equil", expand->elmceq, elmceq_names);
1724     ITYPE ("weight-equil-number-all-lambda", expand->equil_n_at_lam, -1);
1725     ITYPE ("weight-equil-number-samples", expand->equil_samples, -1);
1726     ITYPE ("weight-equil-number-steps", expand->equil_steps, -1);
1727     RTYPE ("weight-equil-wl-delta", expand->equil_wl_delta, -1);
1728     RTYPE ("weight-equil-count-ratio", expand->equil_ratio, -1);
1729     CCTYPE("Seed for Monte Carlo in lambda space");
1730     ITYPE ("lmc-seed", expand->lmc_seed, -1);
1731     RTYPE ("mc-temperature", expand->mc_temp, -1);
1732     ITYPE ("lmc-repeats", expand->lmc_repeats, 1);
1733     ITYPE ("lmc-gibbsdelta", expand->gibbsdeltalam, -1);
1734     ITYPE ("lmc-forced-nstart", expand->lmc_forced_nstart, 0);
1735     EETYPE("symmetrized-transition-matrix", expand->bSymmetrizedTMatrix, yesno_names);
1736     ITYPE("nst-transition-matrix", expand->nstTij, -1);
1737     ITYPE ("mininum-var-min", expand->minvarmin, 100); /*default is reasonable */
1738     ITYPE ("weight-c-range", expand->c_range, 0);      /* default is just C=0 */
1739     RTYPE ("wl-scale", expand->wl_scale, 0.8);
1740     RTYPE ("wl-ratio", expand->wl_ratio, 0.8);
1741     RTYPE ("init-wl-delta", expand->init_wl_delta, 1.0);
1742     EETYPE("wl-oneovert", expand->bWLoneovert, yesno_names);
1743
1744     *ninp_p   = ninp;
1745     *inp_p    = inp;
1746
1747     return;
1748 }
1749
1750 void get_ir(const char *mdparin, const char *mdparout,
1751             t_inputrec *ir, t_gromppopts *opts,
1752             warninp_t wi)
1753 {
1754     char       *dumstr[2];
1755     double      dumdub[2][6];
1756     t_inpfile  *inp;
1757     const char *tmp;
1758     int         i, j, m, ninp;
1759     char        warn_buf[STRLEN];
1760     t_lambda   *fep    = ir->fepvals;
1761     t_expanded *expand = ir->expandedvals;
1762
1763     init_inputrec_strings();
1764     inp = read_inpfile(mdparin, &ninp, wi);
1765
1766     snew(dumstr[0], STRLEN);
1767     snew(dumstr[1], STRLEN);
1768
1769     if (-1 == search_einp(ninp, inp, "cutoff-scheme"))
1770     {
1771         sprintf(warn_buf,
1772                 "%s did not specify a value for the .mdp option "
1773                 "\"cutoff-scheme\". Probably it was first intended for use "
1774                 "with GROMACS before 4.6. In 4.6, the Verlet scheme was "
1775                 "introduced, but the group scheme was still the default. "
1776                 "The default is now the Verlet scheme, so you will observe "
1777                 "different behaviour.", mdparin);
1778         warning_note(wi, warn_buf);
1779     }
1780
1781     /* ignore the following deprecated commands */
1782     REM_TYPE("title");
1783     REM_TYPE("cpp");
1784     REM_TYPE("domain-decomposition");
1785     REM_TYPE("andersen-seed");
1786     REM_TYPE("dihre");
1787     REM_TYPE("dihre-fc");
1788     REM_TYPE("dihre-tau");
1789     REM_TYPE("nstdihreout");
1790     REM_TYPE("nstcheckpoint");
1791     REM_TYPE("optimize-fft");
1792
1793     /* replace the following commands with the clearer new versions*/
1794     REPL_TYPE("unconstrained-start", "continuation");
1795     REPL_TYPE("foreign-lambda", "fep-lambdas");
1796     REPL_TYPE("verlet-buffer-drift", "verlet-buffer-tolerance");
1797     REPL_TYPE("nstxtcout", "nstxout-compressed");
1798     REPL_TYPE("xtc-grps", "compressed-x-grps");
1799     REPL_TYPE("xtc-precision", "compressed-x-precision");
1800
1801     CCTYPE ("VARIOUS PREPROCESSING OPTIONS");
1802     CTYPE ("Preprocessor information: use cpp syntax.");
1803     CTYPE ("e.g.: -I/home/joe/doe -I/home/mary/roe");
1804     STYPE ("include", opts->include,  NULL);
1805     CTYPE ("e.g.: -DPOSRES -DFLEXIBLE (note these variable names are case sensitive)");
1806     STYPE ("define",  opts->define,   NULL);
1807
1808     CCTYPE ("RUN CONTROL PARAMETERS");
1809     EETYPE("integrator",  ir->eI,         ei_names);
1810     CTYPE ("Start time and timestep in ps");
1811     RTYPE ("tinit",   ir->init_t, 0.0);
1812     RTYPE ("dt",      ir->delta_t,    0.001);
1813     STEPTYPE ("nsteps",   ir->nsteps,     0);
1814     CTYPE ("For exact run continuation or redoing part of a run");
1815     STEPTYPE ("init-step", ir->init_step,  0);
1816     CTYPE ("Part index is updated automatically on checkpointing (keeps files separate)");
1817     ITYPE ("simulation-part", ir->simulation_part, 1);
1818     CTYPE ("mode for center of mass motion removal");
1819     EETYPE("comm-mode",   ir->comm_mode,  ecm_names);
1820     CTYPE ("number of steps for center of mass motion removal");
1821     ITYPE ("nstcomm", ir->nstcomm,    100);
1822     CTYPE ("group(s) for center of mass motion removal");
1823     STYPE ("comm-grps",   is->vcm,            NULL);
1824
1825     CCTYPE ("LANGEVIN DYNAMICS OPTIONS");
1826     CTYPE ("Friction coefficient (amu/ps) and random seed");
1827     RTYPE ("bd-fric",     ir->bd_fric,    0.0);
1828     STEPTYPE ("ld-seed",  ir->ld_seed,    -1);
1829
1830     /* Em stuff */
1831     CCTYPE ("ENERGY MINIMIZATION OPTIONS");
1832     CTYPE ("Force tolerance and initial step-size");
1833     RTYPE ("emtol",       ir->em_tol,     10.0);
1834     RTYPE ("emstep",      ir->em_stepsize, 0.01);
1835     CTYPE ("Max number of iterations in relax-shells");
1836     ITYPE ("niter",       ir->niter,      20);
1837     CTYPE ("Step size (ps^2) for minimization of flexible constraints");
1838     RTYPE ("fcstep",      ir->fc_stepsize, 0);
1839     CTYPE ("Frequency of steepest descents steps when doing CG");
1840     ITYPE ("nstcgsteep",  ir->nstcgsteep, 1000);
1841     ITYPE ("nbfgscorr",   ir->nbfgscorr,  10);
1842
1843     CCTYPE ("TEST PARTICLE INSERTION OPTIONS");
1844     RTYPE ("rtpi",    ir->rtpi,   0.05);
1845
1846     /* Output options */
1847     CCTYPE ("OUTPUT CONTROL OPTIONS");
1848     CTYPE ("Output frequency for coords (x), velocities (v) and forces (f)");
1849     ITYPE ("nstxout", ir->nstxout,    0);
1850     ITYPE ("nstvout", ir->nstvout,    0);
1851     ITYPE ("nstfout", ir->nstfout,    0);
1852     CTYPE ("Output frequency for energies to log file and energy file");
1853     ITYPE ("nstlog",  ir->nstlog, 1000);
1854     ITYPE ("nstcalcenergy", ir->nstcalcenergy, 100);
1855     ITYPE ("nstenergy",   ir->nstenergy,  1000);
1856     CTYPE ("Output frequency and precision for .xtc file");
1857     ITYPE ("nstxout-compressed", ir->nstxout_compressed,  0);
1858     RTYPE ("compressed-x-precision", ir->x_compression_precision, 1000.0);
1859     CTYPE ("This selects the subset of atoms for the compressed");
1860     CTYPE ("trajectory file. You can select multiple groups. By");
1861     CTYPE ("default, all atoms will be written.");
1862     STYPE ("compressed-x-grps", is->x_compressed_groups, NULL);
1863     CTYPE ("Selection of energy groups");
1864     STYPE ("energygrps",  is->energy,         NULL);
1865
1866     /* Neighbor searching */
1867     CCTYPE ("NEIGHBORSEARCHING PARAMETERS");
1868     CTYPE ("cut-off scheme (Verlet: particle based cut-offs, group: using charge groups)");
1869     EETYPE("cutoff-scheme",     ir->cutoff_scheme,    ecutscheme_names);
1870     CTYPE ("nblist update frequency");
1871     ITYPE ("nstlist", ir->nstlist,    10);
1872     CTYPE ("ns algorithm (simple or grid)");
1873     EETYPE("ns-type",     ir->ns_type,    ens_names);
1874     CTYPE ("Periodic boundary conditions: xyz, no, xy");
1875     EETYPE("pbc",         ir->ePBC,       epbc_names);
1876     EETYPE("periodic-molecules", ir->bPeriodicMols, yesno_names);
1877     CTYPE ("Allowed energy error due to the Verlet buffer in kJ/mol/ps per atom,");
1878     CTYPE ("a value of -1 means: use rlist");
1879     RTYPE("verlet-buffer-tolerance", ir->verletbuf_tol,    0.005);
1880     CTYPE ("nblist cut-off");
1881     RTYPE ("rlist",   ir->rlist,  1.0);
1882     CTYPE ("long-range cut-off for switched potentials");
1883     RTYPE ("rlistlong",   ir->rlistlong,  -1);
1884     ITYPE ("nstcalclr",   ir->nstcalclr,  -1);
1885
1886     /* Electrostatics */
1887     CCTYPE ("OPTIONS FOR ELECTROSTATICS AND VDW");
1888     CTYPE ("Method for doing electrostatics");
1889     EETYPE("coulombtype", ir->coulombtype,    eel_names);
1890     EETYPE("coulomb-modifier",    ir->coulomb_modifier,    eintmod_names);
1891     CTYPE ("cut-off lengths");
1892     RTYPE ("rcoulomb-switch", ir->rcoulomb_switch,    0.0);
1893     RTYPE ("rcoulomb",    ir->rcoulomb,   1.0);
1894     CTYPE ("Relative dielectric constant for the medium and the reaction field");
1895     RTYPE ("epsilon-r",   ir->epsilon_r,  1.0);
1896     RTYPE ("epsilon-rf",  ir->epsilon_rf, 0.0);
1897     CTYPE ("Method for doing Van der Waals");
1898     EETYPE("vdw-type",    ir->vdwtype,    evdw_names);
1899     EETYPE("vdw-modifier",    ir->vdw_modifier,    eintmod_names);
1900     CTYPE ("cut-off lengths");
1901     RTYPE ("rvdw-switch", ir->rvdw_switch,    0.0);
1902     RTYPE ("rvdw",    ir->rvdw,   1.0);
1903     CTYPE ("Apply long range dispersion corrections for Energy and Pressure");
1904     EETYPE("DispCorr",    ir->eDispCorr,  edispc_names);
1905     CTYPE ("Extension of the potential lookup tables beyond the cut-off");
1906     RTYPE ("table-extension", ir->tabext, 1.0);
1907     CTYPE ("Separate tables between energy group pairs");
1908     STYPE ("energygrp-table", is->egptable,   NULL);
1909     CTYPE ("Spacing for the PME/PPPM FFT grid");
1910     RTYPE ("fourierspacing", ir->fourier_spacing, 0.12);
1911     CTYPE ("FFT grid size, when a value is 0 fourierspacing will be used");
1912     ITYPE ("fourier-nx",  ir->nkx,         0);
1913     ITYPE ("fourier-ny",  ir->nky,         0);
1914     ITYPE ("fourier-nz",  ir->nkz,         0);
1915     CTYPE ("EWALD/PME/PPPM parameters");
1916     ITYPE ("pme-order",   ir->pme_order,   4);
1917     RTYPE ("ewald-rtol",  ir->ewald_rtol, 0.00001);
1918     RTYPE ("ewald-rtol-lj", ir->ewald_rtol_lj, 0.001);
1919     EETYPE("lj-pme-comb-rule", ir->ljpme_combination_rule, eljpme_names);
1920     EETYPE("ewald-geometry", ir->ewald_geometry, eewg_names);
1921     RTYPE ("epsilon-surface", ir->epsilon_surface, 0.0);
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 = -1, 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 }