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