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