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