src/gromacs/gmxpreprocess/readir.c

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