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