src/gromacs/gmxpreprocess/readir.c

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