src/gromacs/gmxpreprocess/readir.cpp

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