src/gromacs/gmxpreprocess/readir.c

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