src/gromacs/gmxpreprocess/readir.cpp

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