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->edHdLPrintEnergy, edHdLPrintEnergy_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_note(wi, "Free energy is turned off, so we will not decouple the molecule listed in your input.");
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 (fep->edHdLPrintEnergy == edHdLPrintEnergyYES)
2368     {
2369         fep->edHdLPrintEnergy = edHdLPrintEnergyTOTAL;
2370         warning_note(wi, "Old option for dhdl-print-energy given: "
2371                      "changing \"yes\" to \"total\"\n");
2372     }
2373
2374     if (ir->bSimTemp && (fep->edHdLPrintEnergy == edHdLPrintEnergyNO))
2375     {
2376         /* always print out the energy to dhdl if we are doing
2377            expanded ensemble, since we need the total energy for
2378            analysis if the temperature is changing. In some
2379            conditions one may only want the potential energy, so
2380            we will allow that if the appropriate mdp setting has
2381            been enabled. Otherwise, total it is:
2382          */
2383         fep->edHdLPrintEnergy = edHdLPrintEnergyTOTAL;
2384     }
2385
2386     if ((ir->efep != efepNO) || ir->bSimTemp)
2387     {
2388         ir->bExpanded = FALSE;
2389         if ((ir->efep == efepEXPANDED) || ir->bSimTemp)
2390         {
2391             ir->bExpanded = TRUE;
2392         }
2393         do_fep_params(ir, is->fep_lambda, is->lambda_weights);
2394         if (ir->bSimTemp) /* done after fep params */
2395         {
2396             do_simtemp_params(ir);
2397         }
2398     }
2399     else
2400     {
2401         ir->fepvals->n_lambda = 0;
2402     }
2403
2404     /* WALL PARAMETERS */
2405
2406     do_wall_params(ir, is->wall_atomtype, is->wall_density, opts);
2407
2408     /* ORIENTATION RESTRAINT PARAMETERS */
2409
2410     if (opts->bOrire && str_nelem(is->orirefitgrp, MAXPTR, NULL) != 1)
2411     {
2412         warning_error(wi, "ERROR: Need one orientation restraint fit group\n");
2413     }
2414
2415     /* DEFORMATION PARAMETERS */
2416
2417     clear_mat(ir->deform);
2418     for (i = 0; i < 6; i++)
2419     {
2420         dumdub[0][i] = 0;
2421     }
2422     m = sscanf(is->deform, "%lf %lf %lf %lf %lf %lf",
2423                &(dumdub[0][0]), &(dumdub[0][1]), &(dumdub[0][2]),
2424                &(dumdub[0][3]), &(dumdub[0][4]), &(dumdub[0][5]));
2425     for (i = 0; i < 3; i++)
2426     {
2427         ir->deform[i][i] = dumdub[0][i];
2428     }
2429     ir->deform[YY][XX] = dumdub[0][3];
2430     ir->deform[ZZ][XX] = dumdub[0][4];
2431     ir->deform[ZZ][YY] = dumdub[0][5];
2432     if (ir->epc != epcNO)
2433     {
2434         for (i = 0; i < 3; i++)
2435         {
2436             for (j = 0; j <= i; j++)
2437             {
2438                 if (ir->deform[i][j] != 0 && ir->compress[i][j] != 0)
2439                 {
2440                     warning_error(wi, "A box element has deform set and compressibility > 0");
2441                 }
2442             }
2443         }
2444         for (i = 0; i < 3; i++)
2445         {
2446             for (j = 0; j < i; j++)
2447             {
2448                 if (ir->deform[i][j] != 0)
2449                 {
2450                     for (m = j; m < DIM; m++)
2451                     {
2452                         if (ir->compress[m][j] != 0)
2453                         {
2454                             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.");
2455                             warning(wi, warn_buf);
2456                         }
2457                     }
2458                 }
2459             }
2460         }
2461     }
2462
2463     /* Ion/water position swapping checks */
2464     if (ir->eSwapCoords != eswapNO)
2465     {
2466         if (ir->swap->nstswap < 1)
2467         {
2468             warning_error(wi, "swap_frequency must be 1 or larger when ion swapping is requested");
2469         }
2470         if (ir->swap->nAverage < 1)
2471         {
2472             warning_error(wi, "coupl_steps must be 1 or larger.\n");
2473         }
2474         if (ir->swap->threshold < 1.0)
2475         {
2476             warning_error(wi, "Ion count threshold must be at least 1.\n");
2477         }
2478     }
2479
2480     sfree(dumstr[0]);
2481     sfree(dumstr[1]);
2482 }
2483
2484 static int search_QMstring(const char *s, int ng, const char *gn[])
2485 {
2486     /* same as normal search_string, but this one searches QM strings */
2487     int i;
2488
2489     for (i = 0; (i < ng); i++)
2490     {
2491         if (gmx_strcasecmp(s, gn[i]) == 0)
2492         {
2493             return i;
2494         }
2495     }
2496
2497     gmx_fatal(FARGS, "this QM method or basisset (%s) is not implemented\n!", s);
2498
2499     return -1;
2500
2501 } /* search_QMstring */
2502
2503 /* We would like gn to be const as well, but C doesn't allow this */
2504 int search_string(const char *s, int ng, char *gn[])
2505 {
2506     int i;
2507
2508     for (i = 0; (i < ng); i++)
2509     {
2510         if (gmx_strcasecmp(s, gn[i]) == 0)
2511         {
2512             return i;
2513         }
2514     }
2515
2516     gmx_fatal(FARGS,
2517               "Group %s referenced in the .mdp file was not found in the index file.\n"
2518               "Group names must match either [moleculetype] names or custom index group\n"
2519               "names, in which case you must supply an index file to the '-n' option\n"
2520               "of grompp.",
2521               s);
2522
2523     return -1;
2524 }
2525
2526 static gmx_bool do_numbering(int natoms, gmx_groups_t *groups, int ng, char *ptrs[],
2527                              t_blocka *block, char *gnames[],
2528                              int gtype, int restnm,
2529                              int grptp, gmx_bool bVerbose,
2530                              warninp_t wi)
2531 {
2532     unsigned short *cbuf;
2533     t_grps         *grps = &(groups->grps[gtype]);
2534     int             i, j, gid, aj, ognr, ntot = 0;
2535     const char     *title;
2536     gmx_bool        bRest;
2537     char            warn_buf[STRLEN];
2538
2539     if (debug)
2540     {
2541         fprintf(debug, "Starting numbering %d groups of type %d\n", ng, gtype);
2542     }
2543
2544     title = gtypes[gtype];
2545
2546     snew(cbuf, natoms);
2547     /* Mark all id's as not set */
2548     for (i = 0; (i < natoms); i++)
2549     {
2550         cbuf[i] = NOGID;
2551     }
2552
2553     snew(grps->nm_ind, ng+1); /* +1 for possible rest group */
2554     for (i = 0; (i < ng); i++)
2555     {
2556         /* Lookup the group name in the block structure */
2557         gid = search_string(ptrs[i], block->nr, gnames);
2558         if ((grptp != egrptpONE) || (i == 0))
2559         {
2560             grps->nm_ind[grps->nr++] = gid;
2561         }
2562         if (debug)
2563         {
2564             fprintf(debug, "Found gid %d for group %s\n", gid, ptrs[i]);
2565         }
2566
2567         /* Now go over the atoms in the group */
2568         for (j = block->index[gid]; (j < block->index[gid+1]); j++)
2569         {
2570
2571             aj = block->a[j];
2572
2573             /* Range checking */
2574             if ((aj < 0) || (aj >= natoms))
2575             {
2576                 gmx_fatal(FARGS, "Invalid atom number %d in indexfile", aj);
2577             }
2578             /* Lookup up the old group number */
2579             ognr = cbuf[aj];
2580             if (ognr != NOGID)
2581             {
2582                 gmx_fatal(FARGS, "Atom %d in multiple %s groups (%d and %d)",
2583                           aj+1, title, ognr+1, i+1);
2584             }
2585             else
2586             {
2587                 /* Store the group number in buffer */
2588                 if (grptp == egrptpONE)
2589                 {
2590                     cbuf[aj] = 0;
2591                 }
2592                 else
2593                 {
2594                     cbuf[aj] = i;
2595                 }
2596                 ntot++;
2597             }
2598         }
2599     }
2600
2601     /* Now check whether we have done all atoms */
2602     bRest = FALSE;
2603     if (ntot != natoms)
2604     {
2605         if (grptp == egrptpALL)
2606         {
2607             gmx_fatal(FARGS, "%d atoms are not part of any of the %s groups",
2608                       natoms-ntot, title);
2609         }
2610         else if (grptp == egrptpPART)
2611         {
2612             sprintf(warn_buf, "%d atoms are not part of any of the %s groups",
2613                     natoms-ntot, title);
2614             warning_note(wi, warn_buf);
2615         }
2616         /* Assign all atoms currently unassigned to a rest group */
2617         for (j = 0; (j < natoms); j++)
2618         {
2619             if (cbuf[j] == NOGID)
2620             {
2621                 cbuf[j] = grps->nr;
2622                 bRest   = TRUE;
2623             }
2624         }
2625         if (grptp != egrptpPART)
2626         {
2627             if (bVerbose)
2628             {
2629                 fprintf(stderr,
2630                         "Making dummy/rest group for %s containing %d elements\n",
2631                         title, natoms-ntot);
2632             }
2633             /* Add group name "rest" */
2634             grps->nm_ind[grps->nr] = restnm;
2635
2636             /* Assign the rest name to all atoms not currently assigned to a group */
2637             for (j = 0; (j < natoms); j++)
2638             {
2639                 if (cbuf[j] == NOGID)
2640                 {
2641                     cbuf[j] = grps->nr;
2642                 }
2643             }
2644             grps->nr++;
2645         }
2646     }
2647
2648     if (grps->nr == 1 && (ntot == 0 || ntot == natoms))
2649     {
2650         /* All atoms are part of one (or no) group, no index required */
2651         groups->ngrpnr[gtype] = 0;
2652         groups->grpnr[gtype]  = NULL;
2653     }
2654     else
2655     {
2656         groups->ngrpnr[gtype] = natoms;
2657         snew(groups->grpnr[gtype], natoms);
2658         for (j = 0; (j < natoms); j++)
2659         {
2660             groups->grpnr[gtype][j] = cbuf[j];
2661         }
2662     }
2663
2664     sfree(cbuf);
2665
2666     return (bRest && grptp == egrptpPART);
2667 }
2668
2669 static void calc_nrdf(gmx_mtop_t *mtop, t_inputrec *ir, char **gnames)
2670 {
2671     t_grpopts              *opts;
2672     gmx_groups_t           *groups;
2673     t_pull                 *pull;
2674     int                     natoms, ai, aj, i, j, d, g, imin, jmin;
2675     t_iatom                *ia;
2676     int                    *nrdf2, *na_vcm, na_tot;
2677     double                 *nrdf_tc, *nrdf_vcm, nrdf_uc, n_sub = 0;
2678     gmx_mtop_atomloop_all_t aloop;
2679     t_atom                 *atom;
2680     int                     mb, mol, ftype, as;
2681     gmx_molblock_t         *molb;
2682     gmx_moltype_t          *molt;
2683
2684     /* Calculate nrdf.
2685      * First calc 3xnr-atoms for each group
2686      * then subtract half a degree of freedom for each constraint
2687      *
2688      * Only atoms and nuclei contribute to the degrees of freedom...
2689      */
2690
2691     opts = &ir->opts;
2692
2693     groups = &mtop->groups;
2694     natoms = mtop->natoms;
2695
2696     /* Allocate one more for a possible rest group */
2697     /* We need to sum degrees of freedom into doubles,
2698      * since floats give too low nrdf's above 3 million atoms.
2699      */
2700     snew(nrdf_tc, groups->grps[egcTC].nr+1);
2701     snew(nrdf_vcm, groups->grps[egcVCM].nr+1);
2702     snew(na_vcm, groups->grps[egcVCM].nr+1);
2703
2704     for (i = 0; i < groups->grps[egcTC].nr; i++)
2705     {
2706         nrdf_tc[i] = 0;
2707     }
2708     for (i = 0; i < groups->grps[egcVCM].nr+1; i++)
2709     {
2710         nrdf_vcm[i] = 0;
2711     }
2712
2713     snew(nrdf2, natoms);
2714     aloop = gmx_mtop_atomloop_all_init(mtop);
2715     while (gmx_mtop_atomloop_all_next(aloop, &i, &atom))
2716     {
2717         nrdf2[i] = 0;
2718         if (atom->ptype == eptAtom || atom->ptype == eptNucleus)
2719         {
2720             g = ggrpnr(groups, egcFREEZE, i);
2721             /* Double count nrdf for particle i */
2722             for (d = 0; d < DIM; d++)
2723             {
2724                 if (opts->nFreeze[g][d] == 0)
2725                 {
2726                     nrdf2[i] += 2;
2727                 }
2728             }
2729             nrdf_tc [ggrpnr(groups, egcTC, i)]  += 0.5*nrdf2[i];
2730             nrdf_vcm[ggrpnr(groups, egcVCM, i)] += 0.5*nrdf2[i];
2731         }
2732     }
2733
2734     as = 0;
2735     for (mb = 0; mb < mtop->nmolblock; mb++)
2736     {
2737         molb = &mtop->molblock[mb];
2738         molt = &mtop->moltype[molb->type];
2739         atom = molt->atoms.atom;
2740         for (mol = 0; mol < molb->nmol; mol++)
2741         {
2742             for (ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++)
2743             {
2744                 ia = molt->ilist[ftype].iatoms;
2745                 for (i = 0; i < molt->ilist[ftype].nr; )
2746                 {
2747                     /* Subtract degrees of freedom for the constraints,
2748                      * if the particles still have degrees of freedom left.
2749                      * If one of the particles is a vsite or a shell, then all
2750                      * constraint motion will go there, but since they do not
2751                      * contribute to the constraints the degrees of freedom do not
2752                      * change.
2753                      */
2754                     ai = as + ia[1];
2755                     aj = as + ia[2];
2756                     if (((atom[ia[1]].ptype == eptNucleus) ||
2757                          (atom[ia[1]].ptype == eptAtom)) &&
2758                         ((atom[ia[2]].ptype == eptNucleus) ||
2759                          (atom[ia[2]].ptype == eptAtom)))
2760                     {
2761                         if (nrdf2[ai] > 0)
2762                         {
2763                             jmin = 1;
2764                         }
2765                         else
2766                         {
2767                             jmin = 2;
2768                         }
2769                         if (nrdf2[aj] > 0)
2770                         {
2771                             imin = 1;
2772                         }
2773                         else
2774                         {
2775                             imin = 2;
2776                         }
2777                         imin       = min(imin, nrdf2[ai]);
2778                         jmin       = min(jmin, nrdf2[aj]);
2779                         nrdf2[ai] -= imin;
2780                         nrdf2[aj] -= jmin;
2781                         nrdf_tc [ggrpnr(groups, egcTC, ai)]  -= 0.5*imin;
2782                         nrdf_tc [ggrpnr(groups, egcTC, aj)]  -= 0.5*jmin;
2783                         nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5*imin;
2784                         nrdf_vcm[ggrpnr(groups, egcVCM, aj)] -= 0.5*jmin;
2785                     }
2786                     ia += interaction_function[ftype].nratoms+1;
2787                     i  += interaction_function[ftype].nratoms+1;
2788                 }
2789             }
2790             ia = molt->ilist[F_SETTLE].iatoms;
2791             for (i = 0; i < molt->ilist[F_SETTLE].nr; )
2792             {
2793                 /* Subtract 1 dof from every atom in the SETTLE */
2794                 for (j = 0; j < 3; j++)
2795                 {
2796                     ai         = as + ia[1+j];
2797                     imin       = min(2, nrdf2[ai]);
2798                     nrdf2[ai] -= imin;
2799                     nrdf_tc [ggrpnr(groups, egcTC, ai)]  -= 0.5*imin;
2800                     nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5*imin;
2801                 }
2802                 ia += 4;
2803                 i  += 4;
2804             }
2805             as += molt->atoms.nr;
2806         }
2807     }
2808
2809     if (ir->ePull == epullCONSTRAINT)
2810     {
2811         /* Correct nrdf for the COM constraints.
2812          * We correct using the TC and VCM group of the first atom
2813          * in the reference and pull group. If atoms in one pull group
2814          * belong to different TC or VCM groups it is anyhow difficult
2815          * to determine the optimal nrdf assignment.
2816          */
2817         pull = ir->pull;
2818
2819         for (i = 0; i < pull->ncoord; i++)
2820         {
2821             imin = 1;
2822
2823             for (j = 0; j < 2; j++)
2824             {
2825                 const t_pull_group *pgrp;
2826
2827                 pgrp = &pull->group[pull->coord[i].group[j]];
2828
2829                 if (pgrp->nat > 0)
2830                 {
2831                     /* Subtract 1/2 dof from each group */
2832                     ai = pgrp->ind[0];
2833                     nrdf_tc [ggrpnr(groups, egcTC, ai)]  -= 0.5*imin;
2834                     nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5*imin;
2835                     if (nrdf_tc[ggrpnr(groups, egcTC, ai)] < 0)
2836                     {
2837                         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)]]);
2838                     }
2839                 }
2840                 else
2841                 {
2842                     /* We need to subtract the whole DOF from group j=1 */
2843                     imin += 1;
2844                 }
2845             }
2846         }
2847     }
2848
2849     if (ir->nstcomm != 0)
2850     {
2851         /* Subtract 3 from the number of degrees of freedom in each vcm group
2852          * when com translation is removed and 6 when rotation is removed
2853          * as well.
2854          */
2855         switch (ir->comm_mode)
2856         {
2857             case ecmLINEAR:
2858                 n_sub = ndof_com(ir);
2859                 break;
2860             case ecmANGULAR:
2861                 n_sub = 6;
2862                 break;
2863             default:
2864                 n_sub = 0;
2865                 gmx_incons("Checking comm_mode");
2866         }
2867
2868         for (i = 0; i < groups->grps[egcTC].nr; i++)
2869         {
2870             /* Count the number of atoms of TC group i for every VCM group */
2871             for (j = 0; j < groups->grps[egcVCM].nr+1; j++)
2872             {
2873                 na_vcm[j] = 0;
2874             }
2875             na_tot = 0;
2876             for (ai = 0; ai < natoms; ai++)
2877             {
2878                 if (ggrpnr(groups, egcTC, ai) == i)
2879                 {
2880                     na_vcm[ggrpnr(groups, egcVCM, ai)]++;
2881                     na_tot++;
2882                 }
2883             }
2884             /* Correct for VCM removal according to the fraction of each VCM
2885              * group present in this TC group.
2886              */
2887             nrdf_uc = nrdf_tc[i];
2888             if (debug)
2889             {
2890                 fprintf(debug, "T-group[%d] nrdf_uc = %g, n_sub = %g\n",
2891                         i, nrdf_uc, n_sub);
2892             }
2893             nrdf_tc[i] = 0;
2894             for (j = 0; j < groups->grps[egcVCM].nr+1; j++)
2895             {
2896                 if (nrdf_vcm[j] > n_sub)
2897                 {
2898                     nrdf_tc[i] += nrdf_uc*((double)na_vcm[j]/(double)na_tot)*
2899                         (nrdf_vcm[j] - n_sub)/nrdf_vcm[j];
2900                 }
2901                 if (debug)
2902                 {
2903                     fprintf(debug, "  nrdf_vcm[%d] = %g, nrdf = %g\n",
2904                             j, nrdf_vcm[j], nrdf_tc[i]);
2905                 }
2906             }
2907         }
2908     }
2909     for (i = 0; (i < groups->grps[egcTC].nr); i++)
2910     {
2911         opts->nrdf[i] = nrdf_tc[i];
2912         if (opts->nrdf[i] < 0)
2913         {
2914             opts->nrdf[i] = 0;
2915         }
2916         fprintf(stderr,
2917                 "Number of degrees of freedom in T-Coupling group %s is %.2f\n",
2918                 gnames[groups->grps[egcTC].nm_ind[i]], opts->nrdf[i]);
2919     }
2920
2921     sfree(nrdf2);
2922     sfree(nrdf_tc);
2923     sfree(nrdf_vcm);
2924     sfree(na_vcm);
2925 }
2926
2927 static void decode_cos(char *s, t_cosines *cosine)
2928 {
2929     char   *t;
2930     char    format[STRLEN], f1[STRLEN];
2931     double  a, phi;
2932     int     i;
2933
2934     t = strdup(s);
2935     trim(t);
2936
2937     cosine->n   = 0;
2938     cosine->a   = NULL;
2939     cosine->phi = NULL;
2940     if (strlen(t))
2941     {
2942         sscanf(t, "%d", &(cosine->n));
2943         if (cosine->n <= 0)
2944         {
2945             cosine->n = 0;
2946         }
2947         else
2948         {
2949             snew(cosine->a, cosine->n);
2950             snew(cosine->phi, cosine->n);
2951
2952             sprintf(format, "%%*d");
2953             for (i = 0; (i < cosine->n); i++)
2954             {
2955                 strcpy(f1, format);
2956                 strcat(f1, "%lf%lf");
2957                 if (sscanf(t, f1, &a, &phi) < 2)
2958                 {
2959                     gmx_fatal(FARGS, "Invalid input for electric field shift: '%s'", t);
2960                 }
2961                 cosine->a[i]   = a;
2962                 cosine->phi[i] = phi;
2963                 strcat(format, "%*lf%*lf");
2964             }
2965         }
2966     }
2967     sfree(t);
2968 }
2969
2970 static gmx_bool do_egp_flag(t_inputrec *ir, gmx_groups_t *groups,
2971                             const char *option, const char *val, int flag)
2972 {
2973     /* The maximum number of energy group pairs would be MAXPTR*(MAXPTR+1)/2.
2974      * But since this is much larger than STRLEN, such a line can not be parsed.
2975      * The real maximum is the number of names that fit in a string: STRLEN/2.
2976      */
2977 #define EGP_MAX (STRLEN/2)
2978     int      nelem, i, j, k, nr;
2979     char    *names[EGP_MAX];
2980     char  ***gnames;
2981     gmx_bool bSet;
2982
2983     gnames = groups->grpname;
2984
2985     nelem = str_nelem(val, EGP_MAX, names);
2986     if (nelem % 2 != 0)
2987     {
2988         gmx_fatal(FARGS, "The number of groups for %s is odd", option);
2989     }
2990     nr   = groups->grps[egcENER].nr;
2991     bSet = FALSE;
2992     for (i = 0; i < nelem/2; i++)
2993     {
2994         j = 0;
2995         while ((j < nr) &&
2996                gmx_strcasecmp(names[2*i], *(gnames[groups->grps[egcENER].nm_ind[j]])))
2997         {
2998             j++;
2999         }
3000         if (j == nr)
3001         {
3002             gmx_fatal(FARGS, "%s in %s is not an energy group\n",
3003                       names[2*i], option);
3004         }
3005         k = 0;
3006         while ((k < nr) &&
3007                gmx_strcasecmp(names[2*i+1], *(gnames[groups->grps[egcENER].nm_ind[k]])))
3008         {
3009             k++;
3010         }
3011         if (k == nr)
3012         {
3013             gmx_fatal(FARGS, "%s in %s is not an energy group\n",
3014                       names[2*i+1], option);
3015         }
3016         if ((j < nr) && (k < nr))
3017         {
3018             ir->opts.egp_flags[nr*j+k] |= flag;
3019             ir->opts.egp_flags[nr*k+j] |= flag;
3020             bSet = TRUE;
3021         }
3022     }
3023
3024     return bSet;
3025 }
3026
3027
3028 static void make_swap_groups(
3029         t_swapcoords *swap,
3030         char         *swapgname,
3031         char         *splitg0name,
3032         char         *splitg1name,
3033         char         *solgname,
3034         t_blocka     *grps,
3035         char        **gnames)
3036 {
3037     int   ig = -1, i = 0, j;
3038     char *splitg;
3039
3040
3041     /* Just a quick check here, more thorough checks are in mdrun */
3042     if (strcmp(splitg0name, splitg1name) == 0)
3043     {
3044         gmx_fatal(FARGS, "The split groups can not both be '%s'.", splitg0name);
3045     }
3046
3047     /* First get the swap group index atoms */
3048     ig        = search_string(swapgname, grps->nr, gnames);
3049     swap->nat = grps->index[ig+1] - grps->index[ig];
3050     if (swap->nat > 0)
3051     {
3052         fprintf(stderr, "Swap group '%s' contains %d atoms.\n", swapgname, swap->nat);
3053         snew(swap->ind, swap->nat);
3054         for (i = 0; i < swap->nat; i++)
3055         {
3056             swap->ind[i] = grps->a[grps->index[ig]+i];
3057         }
3058     }
3059     else
3060     {
3061         gmx_fatal(FARGS, "You defined an empty group of atoms for swapping.");
3062     }
3063
3064     /* Now do so for the split groups */
3065     for (j = 0; j < 2; j++)
3066     {
3067         if (j == 0)
3068         {
3069             splitg = splitg0name;
3070         }
3071         else
3072         {
3073             splitg = splitg1name;
3074         }
3075
3076         ig                 = search_string(splitg, grps->nr, gnames);
3077         swap->nat_split[j] = grps->index[ig+1] - grps->index[ig];
3078         if (swap->nat_split[j] > 0)
3079         {
3080             fprintf(stderr, "Split group %d '%s' contains %d atom%s.\n",
3081                     j, splitg, swap->nat_split[j], (swap->nat_split[j] > 1) ? "s" : "");
3082             snew(swap->ind_split[j], swap->nat_split[j]);
3083             for (i = 0; i < swap->nat_split[j]; i++)
3084             {
3085                 swap->ind_split[j][i] = grps->a[grps->index[ig]+i];
3086             }
3087         }
3088         else
3089         {
3090             gmx_fatal(FARGS, "Split group %d has to contain at least 1 atom!", j);
3091         }
3092     }
3093
3094     /* Now get the solvent group index atoms */
3095     ig            = search_string(solgname, grps->nr, gnames);
3096     swap->nat_sol = grps->index[ig+1] - grps->index[ig];
3097     if (swap->nat_sol > 0)
3098     {
3099         fprintf(stderr, "Solvent group '%s' contains %d atoms.\n", solgname, swap->nat_sol);
3100         snew(swap->ind_sol, swap->nat_sol);
3101         for (i = 0; i < swap->nat_sol; i++)
3102         {
3103             swap->ind_sol[i] = grps->a[grps->index[ig]+i];
3104         }
3105     }
3106     else
3107     {
3108         gmx_fatal(FARGS, "You defined an empty group of solvent. Cannot exchange ions.");
3109     }
3110 }
3111
3112
3113 void make_IMD_group(t_IMD *IMDgroup, char *IMDgname, t_blocka *grps, char **gnames)
3114 {
3115     int      ig = -1, i;
3116
3117
3118     ig            = search_string(IMDgname, grps->nr, gnames);
3119     IMDgroup->nat = grps->index[ig+1] - grps->index[ig];
3120
3121     if (IMDgroup->nat > 0)
3122     {
3123         fprintf(stderr, "Group '%s' with %d atoms can be activated for interactive molecular dynamics (IMD).\n",
3124                 IMDgname, IMDgroup->nat);
3125         snew(IMDgroup->ind, IMDgroup->nat);
3126         for (i = 0; i < IMDgroup->nat; i++)
3127         {
3128             IMDgroup->ind[i] = grps->a[grps->index[ig]+i];
3129         }
3130     }
3131 }
3132
3133
3134 void do_index(const char* mdparin, const char *ndx,
3135               gmx_mtop_t *mtop,
3136               gmx_bool bVerbose,
3137               t_inputrec *ir, rvec *v,
3138               warninp_t wi)
3139 {
3140     t_blocka     *grps;
3141     gmx_groups_t *groups;
3142     int           natoms;
3143     t_symtab     *symtab;
3144     t_atoms       atoms_all;
3145     char          warnbuf[STRLEN], **gnames;
3146     int           nr, ntcg, ntau_t, nref_t, nacc, nofg, nSA, nSA_points, nSA_time, nSA_temp;
3147     real          tau_min;
3148     int           nstcmin;
3149     int           nacg, nfreeze, nfrdim, nenergy, nvcm, nuser;
3150     char         *ptr1[MAXPTR], *ptr2[MAXPTR], *ptr3[MAXPTR];
3151     int           i, j, k, restnm;
3152     real          SAtime;
3153     gmx_bool      bExcl, bTable, bSetTCpar, bAnneal, bRest;
3154     int           nQMmethod, nQMbasis, nQMcharge, nQMmult, nbSH, nCASorb, nCASelec,
3155                   nSAon, nSAoff, nSAsteps, nQMg, nbOPT, nbTS;
3156     char          warn_buf[STRLEN];
3157
3158     if (bVerbose)
3159     {
3160         fprintf(stderr, "processing index file...\n");
3161     }
3162     debug_gmx();
3163     if (ndx == NULL)
3164     {
3165         snew(grps, 1);
3166         snew(grps->index, 1);
3167         snew(gnames, 1);
3168         atoms_all = gmx_mtop_global_atoms(mtop);
3169         analyse(&atoms_all, grps, &gnames, FALSE, TRUE);
3170         free_t_atoms(&atoms_all, FALSE);
3171     }
3172     else
3173     {
3174         grps = init_index(ndx, &gnames);
3175     }
3176
3177     groups = &mtop->groups;
3178     natoms = mtop->natoms;
3179     symtab = &mtop->symtab;
3180
3181     snew(groups->grpname, grps->nr+1);
3182
3183     for (i = 0; (i < grps->nr); i++)
3184     {
3185         groups->grpname[i] = put_symtab(symtab, gnames[i]);
3186     }
3187     groups->grpname[i] = put_symtab(symtab, "rest");
3188     restnm             = i;
3189     srenew(gnames, grps->nr+1);
3190     gnames[restnm]   = *(groups->grpname[i]);
3191     groups->ngrpname = grps->nr+1;
3192
3193     set_warning_line(wi, mdparin, -1);
3194
3195     ntau_t = str_nelem(is->tau_t, MAXPTR, ptr1);
3196     nref_t = str_nelem(is->ref_t, MAXPTR, ptr2);
3197     ntcg   = str_nelem(is->tcgrps, MAXPTR, ptr3);
3198     if ((ntau_t != ntcg) || (nref_t != ntcg))
3199     {
3200         gmx_fatal(FARGS, "Invalid T coupling input: %d groups, %d ref-t values and "
3201                   "%d tau-t values", ntcg, nref_t, ntau_t);
3202     }
3203
3204     bSetTCpar = (ir->etc || EI_SD(ir->eI) || ir->eI == eiBD || EI_TPI(ir->eI));
3205     do_numbering(natoms, groups, ntcg, ptr3, grps, gnames, egcTC,
3206                  restnm, bSetTCpar ? egrptpALL : egrptpALL_GENREST, bVerbose, wi);
3207     nr            = groups->grps[egcTC].nr;
3208     ir->opts.ngtc = nr;
3209     snew(ir->opts.nrdf, nr);
3210     snew(ir->opts.tau_t, nr);
3211     snew(ir->opts.ref_t, nr);
3212     if (ir->eI == eiBD && ir->bd_fric == 0)
3213     {
3214         fprintf(stderr, "bd-fric=0, so tau-t will be used as the inverse friction constant(s)\n");
3215     }
3216
3217     if (bSetTCpar)
3218     {
3219         if (nr != nref_t)
3220         {
3221             gmx_fatal(FARGS, "Not enough ref-t and tau-t values!");
3222         }
3223
3224         tau_min = 1e20;
3225         for (i = 0; (i < nr); i++)
3226         {
3227             ir->opts.tau_t[i] = strtod(ptr1[i], NULL);
3228             if ((ir->eI == eiBD || ir->eI == eiSD2) && ir->opts.tau_t[i] <= 0)
3229             {
3230                 sprintf(warn_buf, "With integrator %s tau-t should be larger than 0", ei_names[ir->eI]);
3231                 warning_error(wi, warn_buf);
3232             }
3233
3234             if (ir->etc != etcVRESCALE && ir->opts.tau_t[i] == 0)
3235             {
3236                 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.");
3237             }
3238
3239             if (ir->opts.tau_t[i] >= 0)
3240             {
3241                 tau_min = min(tau_min, ir->opts.tau_t[i]);
3242             }
3243         }
3244         if (ir->etc != etcNO && ir->nsttcouple == -1)
3245         {
3246             ir->nsttcouple = ir_optimal_nsttcouple(ir);
3247         }
3248
3249         if (EI_VV(ir->eI))
3250         {
3251             if ((ir->etc == etcNOSEHOOVER) && (ir->epc == epcBERENDSEN))
3252             {
3253                 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");
3254             }
3255             if ((ir->epc == epcMTTK) && (ir->etc > etcNO))
3256             {
3257                 if (ir->nstpcouple != ir->nsttcouple)
3258                 {
3259                     int mincouple = min(ir->nstpcouple, ir->nsttcouple);
3260                     ir->nstpcouple = ir->nsttcouple = mincouple;
3261                     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);
3262                     warning_note(wi, warn_buf);
3263                 }
3264             }
3265         }
3266         /* velocity verlet with averaged kinetic energy KE = 0.5*(v(t+1/2) - v(t-1/2)) is implemented
3267            primarily for testing purposes, and does not work with temperature coupling other than 1 */
3268
3269         if (ETC_ANDERSEN(ir->etc))
3270         {
3271             if (ir->nsttcouple != 1)
3272             {
3273                 ir->nsttcouple = 1;
3274                 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");
3275                 warning_note(wi, warn_buf);
3276             }
3277         }
3278         nstcmin = tcouple_min_integration_steps(ir->etc);
3279         if (nstcmin > 1)
3280         {
3281             if (tau_min/(ir->delta_t*ir->nsttcouple) < nstcmin)
3282             {
3283                 sprintf(warn_buf, "For proper integration of the %s thermostat, tau-t (%g) should be at least %d times larger than nsttcouple*dt (%g)",
3284                         ETCOUPLTYPE(ir->etc),
3285                         tau_min, nstcmin,
3286                         ir->nsttcouple*ir->delta_t);
3287                 warning(wi, warn_buf);
3288             }
3289         }
3290         for (i = 0; (i < nr); i++)
3291         {
3292             ir->opts.ref_t[i] = strtod(ptr2[i], NULL);
3293             if (ir->opts.ref_t[i] < 0)
3294             {
3295                 gmx_fatal(FARGS, "ref-t for group %d negative", i);
3296             }
3297         }
3298         /* set the lambda mc temperature to the md integrator temperature (which should be defined
3299            if we are in this conditional) if mc_temp is negative */
3300         if (ir->expandedvals->mc_temp < 0)
3301         {
3302             ir->expandedvals->mc_temp = ir->opts.ref_t[0]; /*for now, set to the first reft */
3303         }
3304     }
3305
3306     /* Simulated annealing for each group. There are nr groups */
3307     nSA = str_nelem(is->anneal, MAXPTR, ptr1);
3308     if (nSA == 1 && (ptr1[0][0] == 'n' || ptr1[0][0] == 'N'))
3309     {
3310         nSA = 0;
3311     }
3312     if (nSA > 0 && nSA != nr)
3313     {
3314         gmx_fatal(FARGS, "Not enough annealing values: %d (for %d groups)\n", nSA, nr);
3315     }
3316     else
3317     {
3318         snew(ir->opts.annealing, nr);
3319         snew(ir->opts.anneal_npoints, nr);
3320         snew(ir->opts.anneal_time, nr);
3321         snew(ir->opts.anneal_temp, nr);
3322         for (i = 0; i < nr; i++)
3323         {
3324             ir->opts.annealing[i]      = eannNO;
3325             ir->opts.anneal_npoints[i] = 0;
3326             ir->opts.anneal_time[i]    = NULL;
3327             ir->opts.anneal_temp[i]    = NULL;
3328         }
3329         if (nSA > 0)
3330         {
3331             bAnneal = FALSE;
3332             for (i = 0; i < nr; i++)
3333             {
3334                 if (ptr1[i][0] == 'n' || ptr1[i][0] == 'N')
3335                 {
3336                     ir->opts.annealing[i] = eannNO;
3337                 }
3338                 else if (ptr1[i][0] == 's' || ptr1[i][0] == 'S')
3339                 {
3340                     ir->opts.annealing[i] = eannSINGLE;
3341                     bAnneal               = TRUE;
3342                 }
3343                 else if (ptr1[i][0] == 'p' || ptr1[i][0] == 'P')
3344                 {
3345                     ir->opts.annealing[i] = eannPERIODIC;
3346                     bAnneal               = TRUE;
3347                 }
3348             }
3349             if (bAnneal)
3350             {
3351                 /* Read the other fields too */
3352                 nSA_points = str_nelem(is->anneal_npoints, MAXPTR, ptr1);
3353                 if (nSA_points != nSA)
3354                 {
3355                     gmx_fatal(FARGS, "Found %d annealing-npoints values for %d groups\n", nSA_points, nSA);
3356                 }
3357                 for (k = 0, i = 0; i < nr; i++)
3358                 {
3359                     ir->opts.anneal_npoints[i] = strtol(ptr1[i], NULL, 10);
3360                     if (ir->opts.anneal_npoints[i] == 1)
3361                     {
3362                         gmx_fatal(FARGS, "Please specify at least a start and an end point for annealing\n");
3363                     }
3364                     snew(ir->opts.anneal_time[i], ir->opts.anneal_npoints[i]);
3365                     snew(ir->opts.anneal_temp[i], ir->opts.anneal_npoints[i]);
3366                     k += ir->opts.anneal_npoints[i];
3367                 }
3368
3369                 nSA_time = str_nelem(is->anneal_time, MAXPTR, ptr1);
3370                 if (nSA_time != k)
3371                 {
3372                     gmx_fatal(FARGS, "Found %d annealing-time values, wanter %d\n", nSA_time, k);
3373                 }
3374                 nSA_temp = str_nelem(is->anneal_temp, MAXPTR, ptr2);
3375                 if (nSA_temp != k)
3376                 {
3377                     gmx_fatal(FARGS, "Found %d annealing-temp values, wanted %d\n", nSA_temp, k);
3378                 }
3379
3380                 for (i = 0, k = 0; i < nr; i++)
3381                 {
3382
3383                     for (j = 0; j < ir->opts.anneal_npoints[i]; j++)
3384                     {
3385                         ir->opts.anneal_time[i][j] = strtod(ptr1[k], NULL);
3386                         ir->opts.anneal_temp[i][j] = strtod(ptr2[k], NULL);
3387                         if (j == 0)
3388                         {
3389                             if (ir->opts.anneal_time[i][0] > (ir->init_t+GMX_REAL_EPS))
3390                             {
3391                                 gmx_fatal(FARGS, "First time point for annealing > init_t.\n");
3392                             }
3393                         }
3394                         else
3395                         {
3396                             /* j>0 */
3397                             if (ir->opts.anneal_time[i][j] < ir->opts.anneal_time[i][j-1])
3398                             {
3399                                 gmx_fatal(FARGS, "Annealing timepoints out of order: t=%f comes after t=%f\n",
3400                                           ir->opts.anneal_time[i][j], ir->opts.anneal_time[i][j-1]);
3401                             }
3402                         }
3403                         if (ir->opts.anneal_temp[i][j] < 0)
3404                         {
3405                             gmx_fatal(FARGS, "Found negative temperature in annealing: %f\n", ir->opts.anneal_temp[i][j]);
3406                         }
3407                         k++;
3408                     }
3409                 }
3410                 /* Print out some summary information, to make sure we got it right */
3411                 for (i = 0, k = 0; i < nr; i++)
3412                 {
3413                     if (ir->opts.annealing[i] != eannNO)
3414                     {
3415                         j = groups->grps[egcTC].nm_ind[i];
3416                         fprintf(stderr, "Simulated annealing for group %s: %s, %d timepoints\n",
3417                                 *(groups->grpname[j]), eann_names[ir->opts.annealing[i]],
3418                                 ir->opts.anneal_npoints[i]);
3419                         fprintf(stderr, "Time (ps)   Temperature (K)\n");
3420                         /* All terms except the last one */
3421                         for (j = 0; j < (ir->opts.anneal_npoints[i]-1); j++)
3422                         {
3423                             fprintf(stderr, "%9.1f      %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3424                         }
3425
3426                         /* Finally the last one */
3427                         j = ir->opts.anneal_npoints[i]-1;
3428                         if (ir->opts.annealing[i] == eannSINGLE)
3429                         {
3430                             fprintf(stderr, "%9.1f-     %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3431                         }
3432                         else
3433                         {
3434                             fprintf(stderr, "%9.1f      %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3435                             if (fabs(ir->opts.anneal_temp[i][j]-ir->opts.anneal_temp[i][0]) > GMX_REAL_EPS)
3436                             {
3437                                 warning_note(wi, "There is a temperature jump when your annealing loops back.\n");
3438                             }
3439                         }
3440                     }
3441                 }
3442             }
3443         }
3444     }
3445
3446     if (ir->ePull != epullNO)
3447     {
3448         make_pull_groups(ir->pull, is->pull_grp, grps, gnames);
3449
3450         make_pull_coords(ir->pull);
3451     }
3452
3453     if (ir->bRot)
3454     {
3455         make_rotation_groups(ir->rot, is->rot_grp, grps, gnames);
3456     }
3457
3458     if (ir->eSwapCoords != eswapNO)
3459     {
3460         make_swap_groups(ir->swap, swapgrp, splitgrp0, splitgrp1, solgrp, grps, gnames);
3461     }
3462
3463     /* Make indices for IMD session */
3464     if (ir->bIMD)
3465     {
3466         make_IMD_group(ir->imd, is->imd_grp, grps, gnames);
3467     }
3468
3469     nacc = str_nelem(is->acc, MAXPTR, ptr1);
3470     nacg = str_nelem(is->accgrps, MAXPTR, ptr2);
3471     if (nacg*DIM != nacc)
3472     {
3473         gmx_fatal(FARGS, "Invalid Acceleration input: %d groups and %d acc. values",
3474                   nacg, nacc);
3475     }
3476     do_numbering(natoms, groups, nacg, ptr2, grps, gnames, egcACC,
3477                  restnm, egrptpALL_GENREST, bVerbose, wi);
3478     nr = groups->grps[egcACC].nr;
3479     snew(ir->opts.acc, nr);
3480     ir->opts.ngacc = nr;
3481
3482     for (i = k = 0; (i < nacg); i++)
3483     {
3484         for (j = 0; (j < DIM); j++, k++)
3485         {
3486             ir->opts.acc[i][j] = strtod(ptr1[k], NULL);
3487         }
3488     }
3489     for (; (i < nr); i++)
3490     {
3491         for (j = 0; (j < DIM); j++)
3492         {
3493             ir->opts.acc[i][j] = 0;
3494         }
3495     }
3496
3497     nfrdim  = str_nelem(is->frdim, MAXPTR, ptr1);
3498     nfreeze = str_nelem(is->freeze, MAXPTR, ptr2);
3499     if (nfrdim != DIM*nfreeze)
3500     {
3501         gmx_fatal(FARGS, "Invalid Freezing input: %d groups and %d freeze values",
3502                   nfreeze, nfrdim);
3503     }
3504     do_numbering(natoms, groups, nfreeze, ptr2, grps, gnames, egcFREEZE,
3505                  restnm, egrptpALL_GENREST, bVerbose, wi);
3506     nr             = groups->grps[egcFREEZE].nr;
3507     ir->opts.ngfrz = nr;
3508     snew(ir->opts.nFreeze, nr);
3509     for (i = k = 0; (i < nfreeze); i++)
3510     {
3511         for (j = 0; (j < DIM); j++, k++)
3512         {
3513             ir->opts.nFreeze[i][j] = (gmx_strncasecmp(ptr1[k], "Y", 1) == 0);
3514             if (!ir->opts.nFreeze[i][j])
3515             {
3516                 if (gmx_strncasecmp(ptr1[k], "N", 1) != 0)
3517                 {
3518                     sprintf(warnbuf, "Please use Y(ES) or N(O) for freezedim only "
3519                             "(not %s)", ptr1[k]);
3520                     warning(wi, warn_buf);
3521                 }
3522             }
3523         }
3524     }
3525     for (; (i < nr); i++)
3526     {
3527         for (j = 0; (j < DIM); j++)
3528         {
3529             ir->opts.nFreeze[i][j] = 0;
3530         }
3531     }
3532
3533     nenergy = str_nelem(is->energy, MAXPTR, ptr1);
3534     do_numbering(natoms, groups, nenergy, ptr1, grps, gnames, egcENER,
3535                  restnm, egrptpALL_GENREST, bVerbose, wi);
3536     add_wall_energrps(groups, ir->nwall, symtab);
3537     ir->opts.ngener = groups->grps[egcENER].nr;
3538     nvcm            = str_nelem(is->vcm, MAXPTR, ptr1);
3539     bRest           =
3540         do_numbering(natoms, groups, nvcm, ptr1, grps, gnames, egcVCM,
3541                      restnm, nvcm == 0 ? egrptpALL_GENREST : egrptpPART, bVerbose, wi);
3542     if (bRest)
3543     {
3544         warning(wi, "Some atoms are not part of any center of mass motion removal group.\n"
3545                 "This may lead to artifacts.\n"
3546                 "In most cases one should use one group for the whole system.");
3547     }
3548
3549     /* Now we have filled the freeze struct, so we can calculate NRDF */
3550     calc_nrdf(mtop, ir, gnames);
3551
3552     if (v && NULL)
3553     {
3554         real fac, ntot = 0;
3555
3556         /* Must check per group! */
3557         for (i = 0; (i < ir->opts.ngtc); i++)
3558         {
3559             ntot += ir->opts.nrdf[i];
3560         }
3561         if (ntot != (DIM*natoms))
3562         {
3563             fac = sqrt(ntot/(DIM*natoms));
3564             if (bVerbose)
3565             {
3566                 fprintf(stderr, "Scaling velocities by a factor of %.3f to account for constraints\n"
3567                         "and removal of center of mass motion\n", fac);
3568             }
3569             for (i = 0; (i < natoms); i++)
3570             {
3571                 svmul(fac, v[i], v[i]);
3572             }
3573         }
3574     }
3575
3576     nuser = str_nelem(is->user1, MAXPTR, ptr1);
3577     do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcUser1,
3578                  restnm, egrptpALL_GENREST, bVerbose, wi);
3579     nuser = str_nelem(is->user2, MAXPTR, ptr1);
3580     do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcUser2,
3581                  restnm, egrptpALL_GENREST, bVerbose, wi);
3582     nuser = str_nelem(is->x_compressed_groups, MAXPTR, ptr1);
3583     do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcCompressedX,
3584                  restnm, egrptpONE, bVerbose, wi);
3585     nofg = str_nelem(is->orirefitgrp, MAXPTR, ptr1);
3586     do_numbering(natoms, groups, nofg, ptr1, grps, gnames, egcORFIT,
3587                  restnm, egrptpALL_GENREST, bVerbose, wi);
3588
3589     /* QMMM input processing */
3590     nQMg          = str_nelem(is->QMMM, MAXPTR, ptr1);
3591     nQMmethod     = str_nelem(is->QMmethod, MAXPTR, ptr2);
3592     nQMbasis      = str_nelem(is->QMbasis, MAXPTR, ptr3);
3593     if ((nQMmethod != nQMg) || (nQMbasis != nQMg))
3594     {
3595         gmx_fatal(FARGS, "Invalid QMMM input: %d groups %d basissets"
3596                   " and %d methods\n", nQMg, nQMbasis, nQMmethod);
3597     }
3598     /* group rest, if any, is always MM! */
3599     do_numbering(natoms, groups, nQMg, ptr1, grps, gnames, egcQMMM,
3600                  restnm, egrptpALL_GENREST, bVerbose, wi);
3601     nr            = nQMg; /*atoms->grps[egcQMMM].nr;*/
3602     ir->opts.ngQM = nQMg;
3603     snew(ir->opts.QMmethod, nr);
3604     snew(ir->opts.QMbasis, nr);
3605     for (i = 0; i < nr; i++)
3606     {
3607         /* input consists of strings: RHF CASSCF PM3 .. These need to be
3608          * converted to the corresponding enum in names.c
3609          */
3610         ir->opts.QMmethod[i] = search_QMstring(ptr2[i], eQMmethodNR,
3611                                                eQMmethod_names);
3612         ir->opts.QMbasis[i]  = search_QMstring(ptr3[i], eQMbasisNR,
3613                                                eQMbasis_names);
3614
3615     }
3616     nQMmult   = str_nelem(is->QMmult, MAXPTR, ptr1);
3617     nQMcharge = str_nelem(is->QMcharge, MAXPTR, ptr2);
3618     nbSH      = str_nelem(is->bSH, MAXPTR, ptr3);
3619     snew(ir->opts.QMmult, nr);
3620     snew(ir->opts.QMcharge, nr);
3621     snew(ir->opts.bSH, nr);
3622
3623     for (i = 0; i < nr; i++)
3624     {
3625         ir->opts.QMmult[i]   = strtol(ptr1[i], NULL, 10);
3626         ir->opts.QMcharge[i] = strtol(ptr2[i], NULL, 10);
3627         ir->opts.bSH[i]      = (gmx_strncasecmp(ptr3[i], "Y", 1) == 0);
3628     }
3629
3630     nCASelec  = str_nelem(is->CASelectrons, MAXPTR, ptr1);
3631     nCASorb   = str_nelem(is->CASorbitals, MAXPTR, ptr2);
3632     snew(ir->opts.CASelectrons, nr);
3633     snew(ir->opts.CASorbitals, nr);
3634     for (i = 0; i < nr; i++)
3635     {
3636         ir->opts.CASelectrons[i] = strtol(ptr1[i], NULL, 10);
3637         ir->opts.CASorbitals[i]  = strtol(ptr2[i], NULL, 10);
3638     }
3639     /* special optimization options */
3640
3641     nbOPT = str_nelem(is->bOPT, MAXPTR, ptr1);
3642     nbTS  = str_nelem(is->bTS, MAXPTR, ptr2);
3643     snew(ir->opts.bOPT, nr);
3644     snew(ir->opts.bTS, nr);
3645     for (i = 0; i < nr; i++)
3646     {
3647         ir->opts.bOPT[i] = (gmx_strncasecmp(ptr1[i], "Y", 1) == 0);
3648         ir->opts.bTS[i]  = (gmx_strncasecmp(ptr2[i], "Y", 1) == 0);
3649     }
3650     nSAon     = str_nelem(is->SAon, MAXPTR, ptr1);
3651     nSAoff    = str_nelem(is->SAoff, MAXPTR, ptr2);
3652     nSAsteps  = str_nelem(is->SAsteps, MAXPTR, ptr3);
3653     snew(ir->opts.SAon, nr);
3654     snew(ir->opts.SAoff, nr);
3655     snew(ir->opts.SAsteps, nr);
3656
3657     for (i = 0; i < nr; i++)
3658     {
3659         ir->opts.SAon[i]    = strtod(ptr1[i], NULL);
3660         ir->opts.SAoff[i]   = strtod(ptr2[i], NULL);
3661         ir->opts.SAsteps[i] = strtol(ptr3[i], NULL, 10);
3662     }
3663     /* end of QMMM input */
3664
3665     if (bVerbose)
3666     {
3667         for (i = 0; (i < egcNR); i++)
3668         {
3669             fprintf(stderr, "%-16s has %d element(s):", gtypes[i], groups->grps[i].nr);
3670             for (j = 0; (j < groups->grps[i].nr); j++)
3671             {
3672                 fprintf(stderr, " %s", *(groups->grpname[groups->grps[i].nm_ind[j]]));
3673             }
3674             fprintf(stderr, "\n");
3675         }
3676     }
3677
3678     nr = groups->grps[egcENER].nr;
3679     snew(ir->opts.egp_flags, nr*nr);
3680
3681     bExcl = do_egp_flag(ir, groups, "energygrp-excl", is->egpexcl, EGP_EXCL);
3682     if (bExcl && ir->cutoff_scheme == ecutsVERLET)
3683     {
3684         warning_error(wi, "Energy group exclusions are not (yet) implemented for the Verlet scheme");
3685     }
3686     if (bExcl && EEL_FULL(ir->coulombtype))
3687     {
3688         warning(wi, "Can not exclude the lattice Coulomb energy between energy groups");
3689     }
3690
3691     bTable = do_egp_flag(ir, groups, "energygrp-table", is->egptable, EGP_TABLE);
3692     if (bTable && !(ir->vdwtype == evdwUSER) &&
3693         !(ir->coulombtype == eelUSER) && !(ir->coulombtype == eelPMEUSER) &&
3694         !(ir->coulombtype == eelPMEUSERSWITCH))
3695     {
3696         gmx_fatal(FARGS, "Can only have energy group pair tables in combination with user tables for VdW and/or Coulomb");
3697     }
3698
3699     decode_cos(is->efield_x, &(ir->ex[XX]));
3700     decode_cos(is->efield_xt, &(ir->et[XX]));
3701     decode_cos(is->efield_y, &(ir->ex[YY]));
3702     decode_cos(is->efield_yt, &(ir->et[YY]));
3703     decode_cos(is->efield_z, &(ir->ex[ZZ]));
3704     decode_cos(is->efield_zt, &(ir->et[ZZ]));
3705
3706     if (ir->bAdress)
3707     {
3708         do_adress_index(ir->adress, groups, gnames, &(ir->opts), wi);
3709     }
3710
3711     for (i = 0; (i < grps->nr); i++)
3712     {
3713         sfree(gnames[i]);
3714     }
3715     sfree(gnames);
3716     done_blocka(grps);
3717     sfree(grps);
3718
3719 }
3720
3721
3722
3723 static void check_disre(gmx_mtop_t *mtop)
3724 {
3725     gmx_ffparams_t *ffparams;
3726     t_functype     *functype;
3727     t_iparams      *ip;
3728     int             i, ndouble, ftype;
3729     int             label, old_label;
3730
3731     if (gmx_mtop_ftype_count(mtop, F_DISRES) > 0)
3732     {
3733         ffparams  = &mtop->ffparams;
3734         functype  = ffparams->functype;
3735         ip        = ffparams->iparams;
3736         ndouble   = 0;
3737         old_label = -1;
3738         for (i = 0; i < ffparams->ntypes; i++)
3739         {
3740             ftype = functype[i];
3741             if (ftype == F_DISRES)
3742             {
3743                 label = ip[i].disres.label;
3744                 if (label == old_label)
3745                 {
3746                     fprintf(stderr, "Distance restraint index %d occurs twice\n", label);
3747                     ndouble++;
3748                 }
3749                 old_label = label;
3750             }
3751         }
3752         if (ndouble > 0)
3753         {
3754             gmx_fatal(FARGS, "Found %d double distance restraint indices,\n"
3755                       "probably the parameters for multiple pairs in one restraint "
3756                       "are not identical\n", ndouble);
3757         }
3758     }
3759 }
3760
3761 static gmx_bool absolute_reference(t_inputrec *ir, gmx_mtop_t *sys,
3762                                    gmx_bool posres_only,
3763                                    ivec AbsRef)
3764 {
3765     int                  d, g, i;
3766     gmx_mtop_ilistloop_t iloop;
3767     t_ilist             *ilist;
3768     int                  nmol;
3769     t_iparams           *pr;
3770
3771     clear_ivec(AbsRef);
3772
3773     if (!posres_only)
3774     {
3775         /* Check the COM */
3776         for (d = 0; d < DIM; d++)
3777         {
3778             AbsRef[d] = (d < ndof_com(ir) ? 0 : 1);
3779         }
3780         /* Check for freeze groups */
3781         for (g = 0; g < ir->opts.ngfrz; g++)
3782         {
3783             for (d = 0; d < DIM; d++)
3784             {
3785                 if (ir->opts.nFreeze[g][d] != 0)
3786                 {
3787                     AbsRef[d] = 1;
3788                 }
3789             }
3790         }
3791     }
3792
3793     /* Check for position restraints */
3794     iloop = gmx_mtop_ilistloop_init(sys);
3795     while (gmx_mtop_ilistloop_next(iloop, &ilist, &nmol))
3796     {
3797         if (nmol > 0 &&
3798             (AbsRef[XX] == 0 || AbsRef[YY] == 0 || AbsRef[ZZ] == 0))
3799         {
3800             for (i = 0; i < ilist[F_POSRES].nr; i += 2)
3801             {
3802                 pr = &sys->ffparams.iparams[ilist[F_POSRES].iatoms[i]];
3803                 for (d = 0; d < DIM; d++)
3804                 {
3805                     if (pr->posres.fcA[d] != 0)
3806                     {
3807                         AbsRef[d] = 1;
3808                     }
3809                 }
3810             }
3811             for (i = 0; i < ilist[F_FBPOSRES].nr; i += 2)
3812             {
3813                 /* Check for flat-bottom posres */
3814                 pr = &sys->ffparams.iparams[ilist[F_FBPOSRES].iatoms[i]];
3815                 if (pr->fbposres.k != 0)
3816                 {
3817                     switch (pr->fbposres.geom)
3818                     {
3819                         case efbposresSPHERE:
3820                             AbsRef[XX] = AbsRef[YY] = AbsRef[ZZ] = 1;
3821                             break;
3822                         case efbposresCYLINDER:
3823                             AbsRef[XX] = AbsRef[YY] = 1;
3824                             break;
3825                         case efbposresX: /* d=XX */
3826                         case efbposresY: /* d=YY */
3827                         case efbposresZ: /* d=ZZ */
3828                             d         = pr->fbposres.geom - efbposresX;
3829                             AbsRef[d] = 1;
3830                             break;
3831                         default:
3832                             gmx_fatal(FARGS, " Invalid geometry for flat-bottom position restraint.\n"
3833                                       "Expected nr between 1 and %d. Found %d\n", efbposresNR-1,
3834                                       pr->fbposres.geom);
3835                     }
3836                 }
3837             }
3838         }
3839     }
3840
3841     return (AbsRef[XX] != 0 && AbsRef[YY] != 0 && AbsRef[ZZ] != 0);
3842 }
3843
3844 static void
3845 check_combination_rule_differences(const gmx_mtop_t *mtop, int state,
3846                                    gmx_bool *bC6ParametersWorkWithGeometricRules,
3847                                    gmx_bool *bC6ParametersWorkWithLBRules,
3848                                    gmx_bool *bLBRulesPossible)
3849 {
3850     int           ntypes, tpi, tpj, thisLBdiff, thisgeomdiff;
3851     int          *typecount;
3852     real          tol;
3853     double        geometricdiff, LBdiff;
3854     double        c6i, c6j, c12i, c12j;
3855     double        c6, c6_geometric, c6_LB;
3856     double        sigmai, sigmaj, epsi, epsj;
3857     gmx_bool      bCanDoLBRules, bCanDoGeometricRules;
3858     const char   *ptr;
3859
3860     /* A tolerance of 1e-5 seems reasonable for (possibly hand-typed)
3861      * force-field floating point parameters.
3862      */
3863     tol = 1e-5;
3864     ptr = getenv("GMX_LJCOMB_TOL");
3865     if (ptr != NULL)
3866     {
3867         double dbl;
3868
3869         sscanf(ptr, "%lf", &dbl);
3870         tol = dbl;
3871     }
3872
3873     *bC6ParametersWorkWithLBRules         = TRUE;
3874     *bC6ParametersWorkWithGeometricRules  = TRUE;
3875     bCanDoLBRules                         = TRUE;
3876     bCanDoGeometricRules                  = TRUE;
3877     ntypes                                = mtop->ffparams.atnr;
3878     snew(typecount, ntypes);
3879     gmx_mtop_count_atomtypes(mtop, state, typecount);
3880     geometricdiff           = LBdiff = 0.0;
3881     *bLBRulesPossible       = TRUE;
3882     for (tpi = 0; tpi < ntypes; ++tpi)
3883     {
3884         c6i  = mtop->ffparams.iparams[(ntypes + 1) * tpi].lj.c6;
3885         c12i = mtop->ffparams.iparams[(ntypes + 1) * tpi].lj.c12;
3886         for (tpj = tpi; tpj < ntypes; ++tpj)
3887         {
3888             c6j          = mtop->ffparams.iparams[(ntypes + 1) * tpj].lj.c6;
3889             c12j         = mtop->ffparams.iparams[(ntypes + 1) * tpj].lj.c12;
3890             c6           = mtop->ffparams.iparams[ntypes * tpi + tpj].lj.c6;
3891             c6_geometric = sqrt(c6i * c6j);
3892             if (!gmx_numzero(c6_geometric))
3893             {
3894                 if (!gmx_numzero(c12i) && !gmx_numzero(c12j))
3895                 {
3896                     sigmai   = pow(c12i / c6i, 1.0/6.0);
3897                     sigmaj   = pow(c12j / c6j, 1.0/6.0);
3898                     epsi     = c6i * c6i /(4.0 * c12i);
3899                     epsj     = c6j * c6j /(4.0 * c12j);
3900                     c6_LB    = 4.0 * pow(epsi * epsj, 1.0/2.0) * pow(0.5 * (sigmai + sigmaj), 6);
3901                 }
3902                 else
3903                 {
3904                     *bLBRulesPossible = FALSE;
3905                     c6_LB             = c6_geometric;
3906                 }
3907                 bCanDoLBRules = gmx_within_tol(c6_LB, c6, tol);
3908             }
3909
3910             if (FALSE == bCanDoLBRules)
3911             {
3912                 *bC6ParametersWorkWithLBRules = FALSE;
3913             }
3914
3915             bCanDoGeometricRules = gmx_within_tol(c6_geometric, c6, tol);
3916
3917             if (FALSE == bCanDoGeometricRules)
3918             {
3919                 *bC6ParametersWorkWithGeometricRules = FALSE;
3920             }
3921         }
3922     }
3923     sfree(typecount);
3924 }
3925
3926 static void
3927 check_combination_rules(const t_inputrec *ir, const gmx_mtop_t *mtop,
3928                         warninp_t wi)
3929 {
3930     char     err_buf[256];
3931     gmx_bool bLBRulesPossible, bC6ParametersWorkWithGeometricRules, bC6ParametersWorkWithLBRules;
3932
3933     check_combination_rule_differences(mtop, 0,
3934                                        &bC6ParametersWorkWithGeometricRules,
3935                                        &bC6ParametersWorkWithLBRules,
3936                                        &bLBRulesPossible);
3937     if (ir->ljpme_combination_rule == eljpmeLB)
3938     {
3939         if (FALSE == bC6ParametersWorkWithLBRules || FALSE == bLBRulesPossible)
3940         {
3941             warning(wi, "You are using arithmetic-geometric combination rules "
3942                     "in LJ-PME, but your non-bonded C6 parameters do not "
3943                     "follow these rules.");
3944         }
3945     }
3946     else
3947     {
3948         if (FALSE == bC6ParametersWorkWithGeometricRules)
3949         {
3950             if (ir->eDispCorr != edispcNO)
3951             {
3952                 warning_note(wi, "You are using geometric combination rules in "
3953                              "LJ-PME, but your non-bonded C6 parameters do "
3954                              "not follow these rules. "
3955                              "This will introduce very small errors in the forces and energies in "
3956                              "your simulations. Dispersion correction will correct total energy "
3957                              "and/or pressure for isotropic systems, but not forces or surface tensions.");
3958             }
3959             else
3960             {
3961                 warning_note(wi, "You are using geometric combination rules in "
3962                              "LJ-PME, but your non-bonded C6 parameters do "
3963                              "not follow these rules. "
3964                              "This will introduce very small errors in the forces and energies in "
3965                              "your simulations. If your system is homogeneous, consider using dispersion correction "
3966                              "for the total energy and pressure.");
3967             }
3968         }
3969     }
3970 }
3971
3972 void triple_check(const char *mdparin, t_inputrec *ir, gmx_mtop_t *sys,
3973                   warninp_t wi)
3974 {
3975     char                      err_buf[STRLEN];
3976     int                       i, m, c, nmol, npct;
3977     gmx_bool                  bCharge, bAcc;
3978     real                      gdt_max, *mgrp, mt;
3979     rvec                      acc;
3980     gmx_mtop_atomloop_block_t aloopb;
3981     gmx_mtop_atomloop_all_t   aloop;
3982     t_atom                   *atom;
3983     ivec                      AbsRef;
3984     char                      warn_buf[STRLEN];
3985
3986     set_warning_line(wi, mdparin, -1);
3987
3988     if (ir->cutoff_scheme == ecutsVERLET &&
3989         ir->verletbuf_tol > 0 &&
3990         ir->nstlist > 1 &&
3991         ((EI_MD(ir->eI) || EI_SD(ir->eI)) &&
3992          (ir->etc == etcVRESCALE || ir->etc == etcBERENDSEN)))
3993     {
3994         /* Check if a too small Verlet buffer might potentially
3995          * cause more drift than the thermostat can couple off.
3996          */
3997         /* Temperature error fraction for warning and suggestion */
3998         const real T_error_warn    = 0.002;
3999         const real T_error_suggest = 0.001;
4000         /* For safety: 2 DOF per atom (typical with constraints) */
4001         const real nrdf_at         = 2;
4002         real       T, tau, max_T_error;
4003         int        i;
4004
4005         T   = 0;
4006         tau = 0;
4007         for (i = 0; i < ir->opts.ngtc; i++)
4008         {
4009             T   = max(T, ir->opts.ref_t[i]);
4010             tau = max(tau, ir->opts.tau_t[i]);
4011         }
4012         if (T > 0)
4013         {
4014             /* This is a worst case estimate of the temperature error,
4015              * assuming perfect buffer estimation and no cancelation
4016              * of errors. The factor 0.5 is because energy distributes
4017              * equally over Ekin and Epot.
4018              */
4019             max_T_error = 0.5*tau*ir->verletbuf_tol/(nrdf_at*BOLTZ*T);
4020             if (max_T_error > T_error_warn)
4021             {
4022                 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.",
4023                         ir->verletbuf_tol, T, tau,
4024                         100*max_T_error,
4025                         100*T_error_suggest,
4026                         ir->verletbuf_tol*T_error_suggest/max_T_error);
4027                 warning(wi, warn_buf);
4028             }
4029         }
4030     }
4031
4032     if (ETC_ANDERSEN(ir->etc))
4033     {
4034         int i;
4035
4036         for (i = 0; i < ir->opts.ngtc; i++)
4037         {
4038             sprintf(err_buf, "all tau_t must currently be equal using Andersen temperature control, violated for group %d", i);
4039             CHECK(ir->opts.tau_t[0] != ir->opts.tau_t[i]);
4040             sprintf(err_buf, "all tau_t must be postive using Andersen temperature control, tau_t[%d]=%10.6f",
4041                     i, ir->opts.tau_t[i]);
4042             CHECK(ir->opts.tau_t[i] < 0);
4043         }
4044
4045         for (i = 0; i < ir->opts.ngtc; i++)
4046         {
4047             int nsteps = (int)(ir->opts.tau_t[i]/ir->delta_t);
4048             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);
4049             CHECK((nsteps % ir->nstcomm) && (ir->etc == etcANDERSENMASSIVE));
4050         }
4051     }
4052
4053     if (EI_DYNAMICS(ir->eI) && !EI_SD(ir->eI) && ir->eI != eiBD &&
4054         ir->comm_mode == ecmNO &&
4055         !(absolute_reference(ir, sys, FALSE, AbsRef) || ir->nsteps <= 10) &&
4056         !ETC_ANDERSEN(ir->etc))
4057     {
4058         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");
4059     }
4060
4061     /* Check for pressure coupling with absolute position restraints */
4062     if (ir->epc != epcNO && ir->refcoord_scaling == erscNO)
4063     {
4064         absolute_reference(ir, sys, TRUE, AbsRef);
4065         {
4066             for (m = 0; m < DIM; m++)
4067             {
4068                 if (AbsRef[m] && norm2(ir->compress[m]) > 0)
4069                 {
4070                     warning(wi, "You are using pressure coupling with absolute position restraints, this will give artifacts. Use the refcoord_scaling option.");
4071                     break;
4072                 }
4073             }
4074         }
4075     }
4076
4077     bCharge = FALSE;
4078     aloopb  = gmx_mtop_atomloop_block_init(sys);
4079     while (gmx_mtop_atomloop_block_next(aloopb, &atom, &nmol))
4080     {
4081         if (atom->q != 0 || atom->qB != 0)
4082         {
4083             bCharge = TRUE;
4084         }
4085     }
4086
4087     if (!bCharge)
4088     {
4089         if (EEL_FULL(ir->coulombtype))
4090         {
4091             sprintf(err_buf,
4092                     "You are using full electrostatics treatment %s for a system without charges.\n"
4093                     "This costs a lot of performance for just processing zeros, consider using %s instead.\n",
4094                     EELTYPE(ir->coulombtype), EELTYPE(eelCUT));
4095             warning(wi, err_buf);
4096         }
4097     }
4098     else
4099     {
4100         if (ir->coulombtype == eelCUT && ir->rcoulomb > 0 && !ir->implicit_solvent)
4101         {
4102             sprintf(err_buf,
4103                     "You are using a plain Coulomb cut-off, which might produce artifacts.\n"
4104                     "You might want to consider using %s electrostatics.\n",
4105                     EELTYPE(eelPME));
4106             warning_note(wi, err_buf);
4107         }
4108     }
4109
4110     /* Check if combination rules used in LJ-PME are the same as in the force field */
4111     if (EVDW_PME(ir->vdwtype))
4112     {
4113         check_combination_rules(ir, sys, wi);
4114     }
4115
4116     /* Generalized reaction field */
4117     if (ir->opts.ngtc == 0)
4118     {
4119         sprintf(err_buf, "No temperature coupling while using coulombtype %s",
4120                 eel_names[eelGRF]);
4121         CHECK(ir->coulombtype == eelGRF);
4122     }
4123     else
4124     {
4125         sprintf(err_buf, "When using coulombtype = %s"
4126                 " ref-t for temperature coupling should be > 0",
4127                 eel_names[eelGRF]);
4128         CHECK((ir->coulombtype == eelGRF) && (ir->opts.ref_t[0] <= 0));
4129     }
4130
4131     if (ir->eI == eiSD2)
4132     {
4133         sprintf(warn_buf, "The stochastic dynamics integrator %s is deprecated, since\n"
4134                 "it is slower than integrator %s and is slightly less accurate\n"
4135                 "with constraints. Use the %s integrator.",
4136                 ei_names[ir->eI], ei_names[eiSD1], ei_names[eiSD1]);
4137         warning_note(wi, warn_buf);
4138     }
4139
4140     bAcc = FALSE;
4141     for (i = 0; (i < sys->groups.grps[egcACC].nr); i++)
4142     {
4143         for (m = 0; (m < DIM); m++)
4144         {
4145             if (fabs(ir->opts.acc[i][m]) > 1e-6)
4146             {
4147                 bAcc = TRUE;
4148             }
4149         }
4150     }
4151     if (bAcc)
4152     {
4153         clear_rvec(acc);
4154         snew(mgrp, sys->groups.grps[egcACC].nr);
4155         aloop = gmx_mtop_atomloop_all_init(sys);
4156         while (gmx_mtop_atomloop_all_next(aloop, &i, &atom))
4157         {
4158             mgrp[ggrpnr(&sys->groups, egcACC, i)] += atom->m;
4159         }
4160         mt = 0.0;
4161         for (i = 0; (i < sys->groups.grps[egcACC].nr); i++)
4162         {
4163             for (m = 0; (m < DIM); m++)
4164             {
4165                 acc[m] += ir->opts.acc[i][m]*mgrp[i];
4166             }
4167             mt += mgrp[i];
4168         }
4169         for (m = 0; (m < DIM); m++)
4170         {
4171             if (fabs(acc[m]) > 1e-6)
4172             {
4173                 const char *dim[DIM] = { "X", "Y", "Z" };
4174                 fprintf(stderr,
4175                         "Net Acceleration in %s direction, will %s be corrected\n",
4176                         dim[m], ir->nstcomm != 0 ? "" : "not");
4177                 if (ir->nstcomm != 0 && m < ndof_com(ir))
4178                 {
4179                     acc[m] /= mt;
4180                     for (i = 0; (i < sys->groups.grps[egcACC].nr); i++)
4181                     {
4182                         ir->opts.acc[i][m] -= acc[m];
4183                     }
4184                 }
4185             }
4186         }
4187         sfree(mgrp);
4188     }
4189
4190     if (ir->efep != efepNO && ir->fepvals->sc_alpha != 0 &&
4191         !gmx_within_tol(sys->ffparams.reppow, 12.0, 10*GMX_DOUBLE_EPS))
4192     {
4193         gmx_fatal(FARGS, "Soft-core interactions are only supported with VdW repulsion power 12");
4194     }
4195
4196     if (ir->ePull != epullNO)
4197     {
4198         gmx_bool bPullAbsoluteRef;
4199
4200         bPullAbsoluteRef = FALSE;
4201         for (i = 0; i < ir->pull->ncoord; i++)
4202         {
4203             bPullAbsoluteRef = bPullAbsoluteRef ||
4204                 ir->pull->coord[i].group[0] == 0 ||
4205                 ir->pull->coord[i].group[1] == 0;
4206         }
4207         if (bPullAbsoluteRef)
4208         {
4209             absolute_reference(ir, sys, FALSE, AbsRef);
4210             for (m = 0; m < DIM; m++)
4211             {
4212                 if (ir->pull->dim[m] && !AbsRef[m])
4213                 {
4214                     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.");
4215                     break;
4216                 }
4217             }
4218         }
4219
4220         if (ir->pull->eGeom == epullgDIRPBC)
4221         {
4222             for (i = 0; i < 3; i++)
4223             {
4224                 for (m = 0; m <= i; m++)
4225                 {
4226                     if ((ir->epc != epcNO && ir->compress[i][m] != 0) ||
4227                         ir->deform[i][m] != 0)
4228                     {
4229                         for (c = 0; c < ir->pull->ncoord; c++)
4230                         {
4231                             if (ir->pull->coord[c].vec[m] != 0)
4232                             {
4233                                 gmx_fatal(FARGS, "Can not have dynamic box while using pull geometry '%s' (dim %c)", EPULLGEOM(ir->pull->eGeom), 'x'+m);
4234                             }
4235                         }
4236                     }
4237                 }
4238             }
4239         }
4240     }
4241
4242     check_disre(sys);
4243 }
4244
4245 void double_check(t_inputrec *ir, matrix box, gmx_bool bConstr, warninp_t wi)
4246 {
4247     real        min_size;
4248     gmx_bool    bTWIN;
4249     char        warn_buf[STRLEN];
4250     const char *ptr;
4251
4252     ptr = check_box(ir->ePBC, box);
4253     if (ptr)
4254     {
4255         warning_error(wi, ptr);
4256     }
4257
4258     if (bConstr && ir->eConstrAlg == econtSHAKE)
4259     {
4260         if (ir->shake_tol <= 0.0)
4261         {
4262             sprintf(warn_buf, "ERROR: shake-tol must be > 0 instead of %g\n",
4263                     ir->shake_tol);
4264             warning_error(wi, warn_buf);
4265         }
4266
4267         if (IR_TWINRANGE(*ir) && ir->nstlist > 1)
4268         {
4269             sprintf(warn_buf, "With twin-range cut-off's and SHAKE the virial and the pressure are incorrect.");
4270             if (ir->epc == epcNO)
4271             {
4272                 warning(wi, warn_buf);
4273             }
4274             else
4275             {
4276                 warning_error(wi, warn_buf);
4277             }
4278         }
4279     }
4280
4281     if ( (ir->eConstrAlg == econtLINCS) && bConstr)
4282     {
4283         /* If we have Lincs constraints: */
4284         if (ir->eI == eiMD && ir->etc == etcNO &&
4285             ir->eConstrAlg == econtLINCS && ir->nLincsIter == 1)
4286         {
4287             sprintf(warn_buf, "For energy conservation with LINCS, lincs_iter should be 2 or larger.\n");
4288             warning_note(wi, warn_buf);
4289         }
4290
4291         if ((ir->eI == eiCG || ir->eI == eiLBFGS) && (ir->nProjOrder < 8))
4292         {
4293             sprintf(warn_buf, "For accurate %s with LINCS constraints, lincs-order should be 8 or more.", ei_names[ir->eI]);
4294             warning_note(wi, warn_buf);
4295         }
4296         if (ir->epc == epcMTTK)
4297         {
4298             warning_error(wi, "MTTK not compatible with lincs -- use shake instead.");
4299         }
4300     }
4301
4302     if (bConstr && ir->epc == epcMTTK)
4303     {
4304         warning_note(wi, "MTTK with constraints is deprecated, and will be removed in GROMACS 5.1");
4305     }
4306
4307     if (ir->LincsWarnAngle > 90.0)
4308     {
4309         sprintf(warn_buf, "lincs-warnangle can not be larger than 90 degrees, setting it to 90.\n");
4310         warning(wi, warn_buf);
4311         ir->LincsWarnAngle = 90.0;
4312     }
4313
4314     if (ir->ePBC != epbcNONE)
4315     {
4316         if (ir->nstlist == 0)
4317         {
4318             warning(wi, "With nstlist=0 atoms are only put into the box at step 0, therefore drifting atoms might cause the simulation to crash.");
4319         }
4320         bTWIN = (ir->rlistlong > ir->rlist);
4321         if (ir->ns_type == ensGRID)
4322         {
4323             if (sqr(ir->rlistlong) >= max_cutoff2(ir->ePBC, box))
4324             {
4325                 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",
4326                         bTWIN ? (ir->rcoulomb == ir->rlistlong ? "rcoulomb" : "rvdw") : "rlist");
4327                 warning_error(wi, warn_buf);
4328             }
4329         }
4330         else
4331         {
4332             min_size = min(box[XX][XX], min(box[YY][YY], box[ZZ][ZZ]));
4333             if (2*ir->rlistlong >= min_size)
4334             {
4335                 sprintf(warn_buf, "ERROR: One of the box lengths is smaller than twice the cut-off length. Increase the box size or decrease rlist.");
4336                 warning_error(wi, warn_buf);
4337                 if (TRICLINIC(box))
4338                 {
4339                     fprintf(stderr, "Grid search might allow larger cut-off's than simple search with triclinic boxes.");
4340                 }
4341             }
4342         }
4343     }
4344 }
4345
4346 void check_chargegroup_radii(const gmx_mtop_t *mtop, const t_inputrec *ir,
4347                              rvec *x,
4348                              warninp_t wi)
4349 {
4350     real rvdw1, rvdw2, rcoul1, rcoul2;
4351     char warn_buf[STRLEN];
4352
4353     calc_chargegroup_radii(mtop, x, &rvdw1, &rvdw2, &rcoul1, &rcoul2);
4354
4355     if (rvdw1 > 0)
4356     {
4357         printf("Largest charge group radii for Van der Waals: %5.3f, %5.3f nm\n",
4358                rvdw1, rvdw2);
4359     }
4360     if (rcoul1 > 0)
4361     {
4362         printf("Largest charge group radii for Coulomb:       %5.3f, %5.3f nm\n",
4363                rcoul1, rcoul2);
4364     }
4365
4366     if (ir->rlist > 0)
4367     {
4368         if (rvdw1  + rvdw2  > ir->rlist ||
4369             rcoul1 + rcoul2 > ir->rlist)
4370         {
4371             sprintf(warn_buf,
4372                     "The sum of the two largest charge group radii (%f) "
4373                     "is larger than rlist (%f)\n",
4374                     max(rvdw1+rvdw2, rcoul1+rcoul2), ir->rlist);
4375             warning(wi, warn_buf);
4376         }
4377         else
4378         {
4379             /* Here we do not use the zero at cut-off macro,
4380              * since user defined interactions might purposely
4381              * not be zero at the cut-off.
4382              */
4383             if (ir_vdw_is_zero_at_cutoff(ir) &&
4384                 rvdw1 + rvdw2 > ir->rlistlong - ir->rvdw)
4385             {
4386                 sprintf(warn_buf, "The sum of the two largest charge group "
4387                         "radii (%f) is larger than %s (%f) - rvdw (%f).\n"
4388                         "With exact cut-offs, better performance can be "
4389                         "obtained with cutoff-scheme = %s, because it "
4390                         "does not use charge groups at all.",
4391                         rvdw1+rvdw2,
4392                         ir->rlistlong > ir->rlist ? "rlistlong" : "rlist",
4393                         ir->rlistlong, ir->rvdw,
4394                         ecutscheme_names[ecutsVERLET]);
4395                 if (ir_NVE(ir))
4396                 {
4397                     warning(wi, warn_buf);
4398                 }
4399                 else
4400                 {
4401                     warning_note(wi, warn_buf);
4402                 }
4403             }
4404             if (ir_coulomb_is_zero_at_cutoff(ir) &&
4405                 rcoul1 + rcoul2 > ir->rlistlong - ir->rcoulomb)
4406             {
4407                 sprintf(warn_buf, "The sum of the two largest charge group radii (%f) is larger than %s (%f) - rcoulomb (%f).\n"
4408                         "With exact cut-offs, better performance can be obtained with cutoff-scheme = %s, because it does not use charge groups at all.",
4409                         rcoul1+rcoul2,
4410                         ir->rlistlong > ir->rlist ? "rlistlong" : "rlist",
4411                         ir->rlistlong, ir->rcoulomb,
4412                         ecutscheme_names[ecutsVERLET]);
4413                 if (ir_NVE(ir))
4414                 {
4415                     warning(wi, warn_buf);
4416                 }
4417                 else
4418                 {
4419                     warning_note(wi, warn_buf);
4420                 }
4421             }
4422         }
4423     }
4424 }