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