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