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