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