src/gromacs/gmxpreprocess/readir.c

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