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