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