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