src/gromacs/mdtypes/inputrec.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-2010, The GROMACS development team.
   6  * Copyright (c) 2012,2014,2015,2016,2017 by the GROMACS development team.
   7  * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
   8  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   9  * and including many others, as listed in the AUTHORS file in the
  10  * top-level source directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #include "gmxpre.h"
  39
  40 #include "inputrec.h"
  41
  42 #include <cstdio>
  43 #include <cstdlib>
  44 #include <cstring>
  45
  46 #include <algorithm>
  47 #include <numeric>
  48
  49 #include "gromacs/math/veccompare.h"
  50 #include "gromacs/math/vecdump.h"
  51 #include "gromacs/mdtypes/awh_params.h"
  52 #include "gromacs/mdtypes/md_enums.h"
  53 #include "gromacs/mdtypes/multipletimestepping.h"
  54 #include "gromacs/mdtypes/pull_params.h"
  55 #include "gromacs/pbcutil/pbc.h"
  56 #include "gromacs/utility/compare.h"
  57 #include "gromacs/utility/cstringutil.h"
  58 #include "gromacs/utility/fatalerror.h"
  59 #include "gromacs/utility/keyvaluetree.h"
  60 #include "gromacs/utility/smalloc.h"
  61 #include "gromacs/utility/snprintf.h"
  62 #include "gromacs/utility/strconvert.h"
  63 #include "gromacs/utility/stringutil.h"
  64 #include "gromacs/utility/textwriter.h"
  65 #include "gromacs/utility/txtdump.h"
  66
  67 //! Macro to select a bool name
  68 #define EBOOL(e) gmx::boolToString(e)
  69
  70 /* Default values for nstcalcenergy, used when the are no other restrictions. */
  71 constexpr int c_defaultNstCalcEnergy = 10;
  72
  73 /* The minimum number of integration steps required for reasonably accurate
  74  * integration of first and second order coupling algorithms.
  75  */
  76 const int nstmin_berendsen_tcouple = 5;
  77 const int nstmin_berendsen_pcouple = 10;
  78 const int nstmin_harmonic          = 20;
  79
  80 /* Default values for T- and P- coupling intervals, used when the are no other
  81  * restrictions.
  82  */
  83 constexpr int c_defaultNstTCouple = 10;
  84 constexpr int c_defaultNstPCouple = 10;
  85
  86 t_inputrec::t_inputrec()
  87 {
  88     // TODO When this memset is removed, remove the suppression of
  89     // gcc -Wno-class-memaccess in a CMakeLists.txt file.
  90     std::memset(this, 0, sizeof(*this)); // NOLINT(bugprone-undefined-memory-manipulation)
  91     snew(fepvals, 1);
  92     snew(expandedvals, 1);
  93     snew(simtempvals, 1);
  94 }
  95
  96 t_inputrec::~t_inputrec()
  97 {
  98     done_inputrec(this);
  99 }
 100
 101 int ir_optimal_nstcalcenergy(const t_inputrec* ir)
 102 {
 103     int nst;
 104
 105     if (ir->nstlist > 0)
 106     {
 107         nst = ir->nstlist;
 108     }
 109     else
 110     {
 111         nst = c_defaultNstCalcEnergy;
 112     }
 113
 114     if (ir->useMts)
 115     {
 116         nst = std::lcm(nst, ir->mtsLevels.back().stepFactor);
 117     }
 118
 119     return nst;
 120 }
 121
 122 int tcouple_min_integration_steps(int etc)
 123 {
 124     int n;
 125
 126     switch (etc)
 127     {
 128         case etcNO: n = 0; break;
 129         case etcBERENDSEN:
 130         case etcYES: n = nstmin_berendsen_tcouple; break;
 131         case etcVRESCALE:
 132             /* V-rescale supports instantaneous rescaling */
 133             n = 0;
 134             break;
 135         case etcNOSEHOOVER: n = nstmin_harmonic; break;
 136         case etcANDERSEN:
 137         case etcANDERSENMASSIVE: n = 1; break;
 138         default: gmx_incons("Unknown etc value");
 139     }
 140
 141     return n;
 142 }
 143
 144 int ir_optimal_nsttcouple(const t_inputrec* ir)
 145 {
 146     int  nmin, nwanted, n;
 147     real tau_min;
 148     int  g;
 149
 150     nmin = tcouple_min_integration_steps(ir->etc);
 151
 152     nwanted = c_defaultNstTCouple;
 153
 154     tau_min = 1e20;
 155     if (ir->etc != etcNO)
 156     {
 157         for (g = 0; g < ir->opts.ngtc; g++)
 158         {
 159             if (ir->opts.tau_t[g] > 0)
 160             {
 161                 tau_min = std::min(tau_min, ir->opts.tau_t[g]);
 162             }
 163         }
 164     }
 165
 166     if (nmin == 0 || ir->delta_t * nwanted <= tau_min)
 167     {
 168         n = nwanted;
 169     }
 170     else
 171     {
 172         n = static_cast<int>(tau_min / (ir->delta_t * nmin) + 0.001);
 173         if (n < 1)
 174         {
 175             n = 1;
 176         }
 177         while (nwanted % n != 0)
 178         {
 179             n--;
 180         }
 181     }
 182
 183     return n;
 184 }
 185
 186 int pcouple_min_integration_steps(int epc)
 187 {
 188     int n;
 189
 190     switch (epc)
 191     {
 192         case epcNO: n = 0; break;
 193         case epcBERENDSEN:
 194         case epcCRESCALE:
 195         case epcISOTROPIC: n = nstmin_berendsen_pcouple; break;
 196         case epcPARRINELLORAHMAN:
 197         case epcMTTK: n = nstmin_harmonic; break;
 198         default: gmx_incons("Unknown epc value");
 199     }
 200
 201     return n;
 202 }
 203
 204 int ir_optimal_nstpcouple(const t_inputrec* ir)
 205 {
 206     const int minIntegrationSteps = pcouple_min_integration_steps(ir->epc);
 207
 208     const int nwanted = c_defaultNstPCouple;
 209
 210     // With multiple time-stepping we can only compute the pressure at slowest steps
 211     const int minNstPCouple = (ir->useMts ? ir->mtsLevels.back().stepFactor : 1);
 212
 213     int n;
 214     if (minIntegrationSteps == 0 || ir->delta_t * nwanted <= ir->tau_p)
 215     {
 216         n = nwanted;
 217     }
 218     else
 219     {
 220         n = static_cast<int>(ir->tau_p / (ir->delta_t * minIntegrationSteps) + 0.001);
 221         if (n < minNstPCouple)
 222         {
 223             n = minNstPCouple;
 224         }
 225         // Without MTS we try to make nstpcouple a "nice" number
 226         if (!ir->useMts)
 227         {
 228             while (nwanted % n != 0)
 229             {
 230                 n--;
 231             }
 232         }
 233     }
 234
 235     // With MTS, nstpcouple should be a multiple of the slowest MTS interval
 236     if (ir->useMts)
 237     {
 238         n = n - (n % minNstPCouple);
 239     }
 240
 241     return n;
 242 }
 243
 244 gmx_bool ir_coulomb_switched(const t_inputrec* ir)
 245 {
 246     return (ir->coulombtype == eelSWITCH || ir->coulombtype == eelSHIFT
 247             || ir->coulombtype == eelPMESWITCH || ir->coulombtype == eelPMEUSERSWITCH
 248             || ir->coulomb_modifier == eintmodPOTSWITCH || ir->coulomb_modifier == eintmodFORCESWITCH);
 249 }
 250
 251 gmx_bool ir_coulomb_is_zero_at_cutoff(const t_inputrec* ir)
 252 {
 253     return (ir->cutoff_scheme == ecutsVERLET || ir_coulomb_switched(ir)
 254             || ir->coulomb_modifier != eintmodNONE || ir->coulombtype == eelRF_ZERO);
 255 }
 256
 257 gmx_bool ir_coulomb_might_be_zero_at_cutoff(const t_inputrec* ir)
 258 {
 259     return (ir_coulomb_is_zero_at_cutoff(ir) || ir->coulombtype == eelUSER || ir->coulombtype == eelPMEUSER);
 260 }
 261
 262 gmx_bool ir_vdw_switched(const t_inputrec* ir)
 263 {
 264     return (ir->vdwtype == evdwSWITCH || ir->vdwtype == evdwSHIFT
 265             || ir->vdw_modifier == eintmodPOTSWITCH || ir->vdw_modifier == eintmodFORCESWITCH);
 266 }
 267
 268 gmx_bool ir_vdw_is_zero_at_cutoff(const t_inputrec* ir)
 269 {
 270     return (ir->cutoff_scheme == ecutsVERLET || ir_vdw_switched(ir) || ir->vdw_modifier != eintmodNONE);
 271 }
 272
 273 gmx_bool ir_vdw_might_be_zero_at_cutoff(const t_inputrec* ir)
 274 {
 275     return (ir_vdw_is_zero_at_cutoff(ir) || ir->vdwtype == evdwUSER);
 276 }
 277
 278 static void done_lambdas(t_lambda* fep)
 279 {
 280     if (fep->n_lambda > 0)
 281     {
 282         for (int i = 0; i < efptNR; i++)
 283         {
 284             sfree(fep->all_lambda[i]);
 285         }
 286     }
 287     sfree(fep->all_lambda);
 288 }
 289
 290 static void done_t_rot(t_rot* rot)
 291 {
 292     if (rot == nullptr)
 293     {
 294         return;
 295     }
 296     if (rot->grp != nullptr)
 297     {
 298         for (int i = 0; i < rot->ngrp; i++)
 299         {
 300             sfree(rot->grp[i].ind);
 301             sfree(rot->grp[i].x_ref);
 302         }
 303         sfree(rot->grp);
 304     }
 305     sfree(rot);
 306 }
 307
 308 void done_inputrec(t_inputrec* ir)
 309 {
 310     sfree(ir->opts.nrdf);
 311     sfree(ir->opts.ref_t);
 312     for (int i = 0; i < ir->opts.ngtc; i++)
 313     {
 314         sfree(ir->opts.anneal_time[i]);
 315         sfree(ir->opts.anneal_temp[i]);
 316     }
 317     sfree(ir->opts.annealing);
 318     sfree(ir->opts.anneal_npoints);
 319     sfree(ir->opts.anneal_time);
 320     sfree(ir->opts.anneal_temp);
 321     sfree(ir->opts.tau_t);
 322     sfree(ir->opts.acc);
 323     sfree(ir->opts.nFreeze);
 324     sfree(ir->opts.egp_flags);
 325     done_lambdas(ir->fepvals);
 326     sfree(ir->fepvals);
 327     sfree(ir->expandedvals);
 328     sfree(ir->simtempvals);
 329
 330     done_t_rot(ir->rot);
 331     delete ir->params;
 332 }
 333
 334 static void pr_grp_opts(FILE* out, int indent, const char* title, const t_grpopts* opts, gmx_bool bMDPformat)
 335 {
 336     int i, m, j;
 337
 338     if (!bMDPformat)
 339     {
 340         fprintf(out, "%s:\n", title);
 341     }
 342
 343     pr_indent(out, indent);
 344     fprintf(out, "nrdf%s", bMDPformat ? " = " : ":");
 345     for (i = 0; (i < opts->ngtc); i++)
 346     {
 347         fprintf(out, "  %10g", opts->nrdf[i]);
 348     }
 349     fprintf(out, "\n");
 350
 351     pr_indent(out, indent);
 352     fprintf(out, "ref-t%s", bMDPformat ? " = " : ":");
 353     for (i = 0; (i < opts->ngtc); i++)
 354     {
 355         fprintf(out, "  %10g", opts->ref_t[i]);
 356     }
 357     fprintf(out, "\n");
 358
 359     pr_indent(out, indent);
 360     fprintf(out, "tau-t%s", bMDPformat ? " = " : ":");
 361     for (i = 0; (i < opts->ngtc); i++)
 362     {
 363         fprintf(out, "  %10g", opts->tau_t[i]);
 364     }
 365     fprintf(out, "\n");
 366
 367     /* Pretty-print the simulated annealing info */
 368     fprintf(out, "annealing%s", bMDPformat ? " = " : ":");
 369     for (i = 0; (i < opts->ngtc); i++)
 370     {
 371         fprintf(out, "  %10s", EANNEAL(opts->annealing[i]));
 372     }
 373     fprintf(out, "\n");
 374
 375     fprintf(out, "annealing-npoints%s", bMDPformat ? " = " : ":");
 376     for (i = 0; (i < opts->ngtc); i++)
 377     {
 378         fprintf(out, "  %10d", opts->anneal_npoints[i]);
 379     }
 380     fprintf(out, "\n");
 381
 382     for (i = 0; (i < opts->ngtc); i++)
 383     {
 384         if (opts->anneal_npoints[i] > 0)
 385         {
 386             fprintf(out, "annealing-time [%d]:\t", i);
 387             for (j = 0; (j < opts->anneal_npoints[i]); j++)
 388             {
 389                 fprintf(out, "  %10.1f", opts->anneal_time[i][j]);
 390             }
 391             fprintf(out, "\n");
 392             fprintf(out, "annealing-temp [%d]:\t", i);
 393             for (j = 0; (j < opts->anneal_npoints[i]); j++)
 394             {
 395                 fprintf(out, "  %10.1f", opts->anneal_temp[i][j]);
 396             }
 397             fprintf(out, "\n");
 398         }
 399     }
 400
 401     pr_indent(out, indent);
 402     fprintf(out, "acc:\t");
 403     for (i = 0; (i < opts->ngacc); i++)
 404     {
 405         for (m = 0; (m < DIM); m++)
 406         {
 407             fprintf(out, "  %10g", opts->acc[i][m]);
 408         }
 409     }
 410     fprintf(out, "\n");
 411
 412     pr_indent(out, indent);
 413     fprintf(out, "nfreeze:");
 414     for (i = 0; (i < opts->ngfrz); i++)
 415     {
 416         for (m = 0; (m < DIM); m++)
 417         {
 418             fprintf(out, "  %10s", opts->nFreeze[i][m] ? "Y" : "N");
 419         }
 420     }
 421     fprintf(out, "\n");
 422
 423
 424     for (i = 0; (i < opts->ngener); i++)
 425     {
 426         pr_indent(out, indent);
 427         fprintf(out, "energygrp-flags[%3d]:", i);
 428         for (m = 0; (m < opts->ngener); m++)
 429         {
 430             fprintf(out, " %d", opts->egp_flags[opts->ngener * i + m]);
 431         }
 432         fprintf(out, "\n");
 433     }
 434
 435     fflush(out);
 436 }
 437
 438 static void pr_matrix(FILE* fp, int indent, const char* title, const rvec* m, gmx_bool bMDPformat)
 439 {
 440     if (bMDPformat)
 441     {
 442         fprintf(fp,
 443                 "%-10s    = %g %g %g %g %g %g\n",
 444                 title,
 445                 m[XX][XX],
 446                 m[YY][YY],
 447                 m[ZZ][ZZ],
 448                 m[XX][YY],
 449                 m[XX][ZZ],
 450                 m[YY][ZZ]);
 451     }
 452     else
 453     {
 454         pr_rvecs(fp, indent, title, m, DIM);
 455     }
 456 }
 457
 458 #define PS(t, s) pr_str(fp, indent, t, s)
 459 #define PI(t, s) pr_int(fp, indent, t, s)
 460 #define PSTEP(t, s) pr_int64(fp, indent, t, s)
 461 #define PR(t, s) pr_real(fp, indent, t, s)
 462 #define PD(t, s) pr_double(fp, indent, t, s)
 463
 464 static void pr_pull_group(FILE* fp, int indent, int g, const t_pull_group* pgrp)
 465 {
 466     pr_indent(fp, indent);
 467     fprintf(fp, "pull-group %d:\n", g);
 468     indent += 2;
 469     pr_ivec_block(fp, indent, "atom", pgrp->ind.data(), pgrp->ind.size(), TRUE);
 470     pr_rvec(fp, indent, "weight", pgrp->weight.data(), pgrp->weight.size(), TRUE);
 471     PI("pbcatom", pgrp->pbcatom);
 472 }
 473
 474 static void pr_pull_coord(FILE* fp, int indent, int c, const t_pull_coord* pcrd)
 475 {
 476     int g;
 477
 478     pr_indent(fp, indent);
 479     fprintf(fp, "pull-coord %d:\n", c);
 480     PS("type", EPULLTYPE(pcrd->eType));
 481     if (pcrd->eType == epullEXTERNAL)
 482     {
 483         PS("potential-provider", pcrd->externalPotentialProvider.c_str());
 484     }
 485     PS("geometry", EPULLGEOM(pcrd->eGeom));
 486     for (g = 0; g < pcrd->ngroup; g++)
 487     {
 488         char buf[STRLEN];
 489
 490         sprintf(buf, "group[%d]", g);
 491         PI(buf, pcrd->group[g]);
 492     }
 493     pr_ivec(fp, indent, "dim", pcrd->dim, DIM, TRUE);
 494     pr_rvec(fp, indent, "origin", pcrd->origin, DIM, TRUE);
 495     pr_rvec(fp, indent, "vec", pcrd->vec, DIM, TRUE);
 496     PS("start", EBOOL(pcrd->bStart));
 497     PR("init", pcrd->init);
 498     PR("rate", pcrd->rate);
 499     PR("k", pcrd->k);
 500     PR("kB", pcrd->kB);
 501 }
 502
 503 static void pr_simtempvals(FILE* fp, int indent, const t_simtemp* simtemp, int n_lambda)
 504 {
 505     PS("simulated-tempering-scaling", ESIMTEMP(simtemp->eSimTempScale));
 506     PR("sim-temp-low", simtemp->simtemp_low);
 507     PR("sim-temp-high", simtemp->simtemp_high);
 508     pr_rvec(fp, indent, "simulated tempering temperatures", simtemp->temperatures, n_lambda, TRUE);
 509 }
 510
 511 static void pr_expandedvals(FILE* fp, int indent, const t_expanded* expand, int n_lambda)
 512 {
 513
 514     PI("nstexpanded", expand->nstexpanded);
 515     PS("lmc-stats", elamstats_names[expand->elamstats]);
 516     PS("lmc-move", elmcmove_names[expand->elmcmove]);
 517     PS("lmc-weights-equil", elmceq_names[expand->elmceq]);
 518     if (expand->elmceq == elmceqNUMATLAM)
 519     {
 520         PI("weight-equil-number-all-lambda", expand->equil_n_at_lam);
 521     }
 522     if (expand->elmceq == elmceqSAMPLES)
 523     {
 524         PI("weight-equil-number-samples", expand->equil_samples);
 525     }
 526     if (expand->elmceq == elmceqSTEPS)
 527     {
 528         PI("weight-equil-number-steps", expand->equil_steps);
 529     }
 530     if (expand->elmceq == elmceqWLDELTA)
 531     {
 532         PR("weight-equil-wl-delta", expand->equil_wl_delta);
 533     }
 534     if (expand->elmceq == elmceqRATIO)
 535     {
 536         PR("weight-equil-count-ratio", expand->equil_ratio);
 537     }
 538     PI("lmc-seed", expand->lmc_seed);
 539     PR("mc-temperature", expand->mc_temp);
 540     PI("lmc-repeats", expand->lmc_repeats);
 541     PI("lmc-gibbsdelta", expand->gibbsdeltalam);
 542     PI("lmc-forced-nstart", expand->lmc_forced_nstart);
 543     PS("symmetrized-transition-matrix", EBOOL(expand->bSymmetrizedTMatrix));
 544     PI("nst-transition-matrix", expand->nstTij);
 545     PI("mininum-var-min", expand->minvarmin); /*default is reasonable */
 546     PI("weight-c-range", expand->c_range);    /* default is just C=0 */
 547     PR("wl-scale", expand->wl_scale);
 548     PR("wl-ratio", expand->wl_ratio);
 549     PR("init-wl-delta", expand->init_wl_delta);
 550     PS("wl-oneovert", EBOOL(expand->bWLoneovert));
 551
 552     pr_indent(fp, indent);
 553     pr_rvec(fp, indent, "init-lambda-weights", expand->init_lambda_weights, n_lambda, TRUE);
 554     PS("init-weights", EBOOL(expand->bInit_weights));
 555 }
 556
 557 static void pr_fepvals(FILE* fp, int indent, const t_lambda* fep, gmx_bool bMDPformat)
 558 {
 559     int i, j;
 560
 561     PR("init-lambda", fep->init_lambda);
 562     PI("init-lambda-state", fep->init_fep_state);
 563     PR("delta-lambda", fep->delta_lambda);
 564     PI("nstdhdl", fep->nstdhdl);
 565
 566     if (!bMDPformat)
 567     {
 568         PI("n-lambdas", fep->n_lambda);
 569     }
 570     if (fep->n_lambda > 0)
 571     {
 572         pr_indent(fp, indent);
 573         fprintf(fp, "separate-dvdl%s\n", bMDPformat ? " = " : ":");
 574         for (i = 0; i < efptNR; i++)
 575         {
 576             fprintf(fp, "%18s = ", efpt_names[i]);
 577             if (fep->separate_dvdl[i])
 578             {
 579                 fprintf(fp, "  TRUE");
 580             }
 581             else
 582             {
 583                 fprintf(fp, "  FALSE");
 584             }
 585             fprintf(fp, "\n");
 586         }
 587         fprintf(fp, "all-lambdas%s\n", bMDPformat ? " = " : ":");
 588         for (i = 0; i < efptNR; i++)
 589         {
 590             fprintf(fp, "%18s = ", efpt_names[i]);
 591             for (j = 0; j < fep->n_lambda; j++)
 592             {
 593                 fprintf(fp, "  %10g", fep->all_lambda[i][j]);
 594             }
 595             fprintf(fp, "\n");
 596         }
 597     }
 598     PI("calc-lambda-neighbors", fep->lambda_neighbors);
 599     PS("dhdl-print-energy", edHdLPrintEnergy_names[fep->edHdLPrintEnergy]);
 600     PR("sc-alpha", fep->sc_alpha);
 601     PI("sc-power", fep->sc_power);
 602     PR("sc-r-power", fep->sc_r_power);
 603     PR("sc-sigma", fep->sc_sigma);
 604     PR("sc-sigma-min", fep->sc_sigma_min);
 605     PS("sc-coul", EBOOL(fep->bScCoul));
 606     PI("dh-hist-size", fep->dh_hist_size);
 607     PD("dh-hist-spacing", fep->dh_hist_spacing);
 608     PS("separate-dhdl-file", SEPDHDLFILETYPE(fep->separate_dhdl_file));
 609     PS("dhdl-derivatives", DHDLDERIVATIVESTYPE(fep->dhdl_derivatives));
 610 };
 611
 612 static void pr_pull(FILE* fp, int indent, const pull_params_t& pull)
 613 {
 614     int g;
 615
 616     PR("pull-cylinder-r", pull.cylinder_r);
 617     PR("pull-constr-tol", pull.constr_tol);
 618     PS("pull-print-COM", EBOOL(pull.bPrintCOM));
 619     PS("pull-print-ref-value", EBOOL(pull.bPrintRefValue));
 620     PS("pull-print-components", EBOOL(pull.bPrintComp));
 621     PI("pull-nstxout", pull.nstxout);
 622     PI("pull-nstfout", pull.nstfout);
 623     PS("pull-pbc-ref-prev-step-com", EBOOL(pull.bSetPbcRefToPrevStepCOM));
 624     PS("pull-xout-average", EBOOL(pull.bXOutAverage));
 625     PS("pull-fout-average", EBOOL(pull.bFOutAverage));
 626     PI("pull-ngroups", pull.ngroup);
 627     for (g = 0; g < pull.ngroup; g++)
 628     {
 629         pr_pull_group(fp, indent, g, &pull.group[g]);
 630     }
 631     PI("pull-ncoords", pull.ncoord);
 632     for (g = 0; g < pull.ncoord; g++)
 633     {
 634         pr_pull_coord(fp, indent, g, &pull.coord[g]);
 635     }
 636 }
 637
 638 static void pr_awh_bias_dim(FILE* fp, int indent, gmx::AwhDimParams* awhDimParams, const char* prefix)
 639 {
 640     pr_indent(fp, indent);
 641     indent++;
 642     fprintf(fp, "%s:\n", prefix);
 643     PS("coord-provider", EAWHCOORDPROVIDER(awhDimParams->eCoordProvider));
 644     PI("coord-index", awhDimParams->coordIndex + 1);
 645     PR("start", awhDimParams->origin);
 646     PR("end", awhDimParams->end);
 647     PR("period", awhDimParams->period);
 648     PR("force-constant", awhDimParams->forceConstant);
 649     PR("diffusion", awhDimParams->diffusion);
 650     PR("cover-diameter", awhDimParams->coverDiameter);
 651 }
 652
 653 static void pr_awh_bias(FILE* fp, int indent, gmx::AwhBiasParams* awhBiasParams, const char* prefix)
 654 {
 655     char opt[STRLEN];
 656
 657     sprintf(opt, "%s-error-init", prefix);
 658     PR(opt, awhBiasParams->errorInitial);
 659     sprintf(opt, "%s-growth", prefix);
 660     PS(opt, EAWHGROWTH(awhBiasParams->eGrowth));
 661     sprintf(opt, "%s-target", prefix);
 662     PS(opt, EAWHTARGET(awhBiasParams->eTarget));
 663     sprintf(opt, "%s-target-beta-scalng", prefix);
 664     PR(opt, awhBiasParams->targetBetaScaling);
 665     sprintf(opt, "%s-target-cutoff", prefix);
 666     PR(opt, awhBiasParams->targetCutoff);
 667     sprintf(opt, "%s-user-data", prefix);
 668     PS(opt, EBOOL(awhBiasParams->bUserData));
 669     sprintf(opt, "%s-share-group", prefix);
 670     PI(opt, awhBiasParams->shareGroup);
 671     sprintf(opt, "%s-equilibrate-histogram", prefix);
 672     PS(opt, EBOOL(awhBiasParams->equilibrateHistogram));
 673     sprintf(opt, "%s-ndim", prefix);
 674     PI(opt, awhBiasParams->ndim);
 675
 676     for (int d = 0; d < awhBiasParams->ndim; d++)
 677     {
 678         char prefixdim[STRLEN];
 679         sprintf(prefixdim, "%s-dim%d", prefix, d + 1);
 680         pr_awh_bias_dim(fp, indent, &awhBiasParams->dimParams[d], prefixdim);
 681     }
 682 }
 683
 684 static void pr_awh(FILE* fp, int indent, gmx::AwhParams* awhParams)
 685 {
 686     PS("awh-potential", EAWHPOTENTIAL(awhParams->ePotential));
 687     PI("awh-seed", awhParams->seed);
 688     PI("awh-nstout", awhParams->nstOut);
 689     PI("awh-nstsample", awhParams->nstSampleCoord);
 690     PI("awh-nsamples-update", awhParams->numSamplesUpdateFreeEnergy);
 691     PS("awh-share-bias-multisim", EBOOL(awhParams->shareBiasMultisim));
 692     PI("awh-nbias", awhParams->numBias);
 693
 694     for (int k = 0; k < awhParams->numBias; k++)
 695     {
 696         auto prefix = gmx::formatString("awh%d", k + 1);
 697         pr_awh_bias(fp, indent, &awhParams->awhBiasParams[k], prefix.c_str());
 698     }
 699 }
 700
 701 static void pr_rotgrp(FILE* fp, int indent, int g, const t_rotgrp* rotg)
 702 {
 703     pr_indent(fp, indent);
 704     fprintf(fp, "rot-group %d:\n", g);
 705     indent += 2;
 706     PS("rot-type", EROTGEOM(rotg->eType));
 707     PS("rot-massw", EBOOL(rotg->bMassW));
 708     pr_ivec_block(fp, indent, "atom", rotg->ind, rotg->nat, TRUE);
 709     pr_rvecs(fp, indent, "x-ref", rotg->x_ref, rotg->nat);
 710     pr_rvec(fp, indent, "rot-vec", rotg->inputVec, DIM, TRUE);
 711     pr_rvec(fp, indent, "rot-pivot", rotg->pivot, DIM, TRUE);
 712     PR("rot-rate", rotg->rate);
 713     PR("rot-k", rotg->k);
 714     PR("rot-slab-dist", rotg->slab_dist);
 715     PR("rot-min-gauss", rotg->min_gaussian);
 716     PR("rot-eps", rotg->eps);
 717     PS("rot-fit-method", EROTFIT(rotg->eFittype));
 718     PI("rot-potfit-nstep", rotg->PotAngle_nstep);
 719     PR("rot-potfit-step", rotg->PotAngle_step);
 720 }
 721
 722 static void pr_rot(FILE* fp, int indent, const t_rot* rot)
 723 {
 724     int g;
 725
 726     PI("rot-nstrout", rot->nstrout);
 727     PI("rot-nstsout", rot->nstsout);
 728     PI("rot-ngroups", rot->ngrp);
 729     for (g = 0; g < rot->ngrp; g++)
 730     {
 731         pr_rotgrp(fp, indent, g, &rot->grp[g]);
 732     }
 733 }
 734
 735
 736 static void pr_swap(FILE* fp, int indent, const t_swapcoords* swap)
 737 {
 738     char str[STRLEN];
 739
 740     /* Enums for better readability of the code */
 741     enum
 742     {
 743         eCompA = 0,
 744         eCompB
 745     };
 746
 747
 748     PI("swap-frequency", swap->nstswap);
 749
 750     /* The split groups that define the compartments */
 751     for (int j = 0; j < 2; j++)
 752     {
 753         snprintf(str, STRLEN, "massw_split%d", j);
 754         PS(str, EBOOL(swap->massw_split[j]));
 755         snprintf(str, STRLEN, "split atoms group %d", j);
 756         pr_ivec_block(fp, indent, str, swap->grp[j].ind, swap->grp[j].nat, TRUE);
 757     }
 758
 759     /* The solvent group */
 760     snprintf(str, STRLEN, "solvent group %s", swap->grp[eGrpSolvent].molname);
 761     pr_ivec_block(fp, indent, str, swap->grp[eGrpSolvent].ind, swap->grp[eGrpSolvent].nat, TRUE);
 762
 763     /* Now print the indices for all the ion groups: */
 764     for (int ig = eSwapFixedGrpNR; ig < swap->ngrp; ig++)
 765     {
 766         snprintf(str, STRLEN, "ion group %s", swap->grp[ig].molname);
 767         pr_ivec_block(fp, indent, str, swap->grp[ig].ind, swap->grp[ig].nat, TRUE);
 768     }
 769
 770     PR("cyl0-r", swap->cyl0r);
 771     PR("cyl0-up", swap->cyl0u);
 772     PR("cyl0-down", swap->cyl0l);
 773     PR("cyl1-r", swap->cyl1r);
 774     PR("cyl1-up", swap->cyl1u);
 775     PR("cyl1-down", swap->cyl1l);
 776     PI("coupl-steps", swap->nAverage);
 777
 778     /* Print the requested ion counts for both compartments */
 779     for (int ic = eCompA; ic <= eCompB; ic++)
 780     {
 781         for (int ig = eSwapFixedGrpNR; ig < swap->ngrp; ig++)
 782         {
 783             snprintf(str, STRLEN, "%s-in-%c", swap->grp[ig].molname, 'A' + ic);
 784             PI(str, swap->grp[ig].nmolReq[ic]);
 785         }
 786     }
 787
 788     PR("threshold", swap->threshold);
 789     PR("bulk-offsetA", swap->bulkOffset[eCompA]);
 790     PR("bulk-offsetB", swap->bulkOffset[eCompB]);
 791 }
 792
 793
 794 static void pr_imd(FILE* fp, int indent, const t_IMD* imd)
 795 {
 796     PI("IMD-atoms", imd->nat);
 797     pr_ivec_block(fp, indent, "atom", imd->ind, imd->nat, TRUE);
 798 }
 799
 800
 801 void pr_inputrec(FILE* fp, int indent, const char* title, const t_inputrec* ir, gmx_bool bMDPformat)
 802 {
 803     const char* infbuf = "inf";
 804
 805     if (available(fp, ir, indent, title))
 806     {
 807         if (!bMDPformat)
 808         {
 809             indent = pr_title(fp, indent, title);
 810         }
 811         /* Try to make this list appear in the same order as the
 812          * options are written in the default mdout.mdp, and with
 813          * the same user-exposed names to facilitate debugging.
 814          */
 815         PS("integrator", EI(ir->eI));
 816         PR("tinit", ir->init_t);
 817         PR("dt", ir->delta_t);
 818         PSTEP("nsteps", ir->nsteps);
 819         PSTEP("init-step", ir->init_step);
 820         PI("simulation-part", ir->simulation_part);
 821         PS("mts", EBOOL(ir->useMts));
 822         if (ir->useMts)
 823         {
 824             for (int mtsIndex = 1; mtsIndex < static_cast<int>(ir->mtsLevels.size()); mtsIndex++)
 825             {
 826                 const auto&       mtsLevel = ir->mtsLevels[mtsIndex];
 827                 const std::string forceKey = gmx::formatString("mts-level%d-forces", mtsIndex + 1);
 828                 std::string       forceGroups;
 829                 for (int i = 0; i < static_cast<int>(gmx::MtsForceGroups::Count); i++)
 830                 {
 831                     if (mtsLevel.forceGroups[i])
 832                     {
 833                         if (!forceGroups.empty())
 834                         {
 835                             forceGroups += " ";
 836                         }
 837                         forceGroups += gmx::mtsForceGroupNames[i];
 838                     }
 839                 }
 840                 PS(forceKey.c_str(), forceGroups.c_str());
 841                 const std::string factorKey = gmx::formatString("mts-level%d-factor", mtsIndex + 1);
 842                 PI(factorKey.c_str(), mtsLevel.stepFactor);
 843             }
 844         }
 845         PS("comm-mode", ECOM(ir->comm_mode));
 846         PI("nstcomm", ir->nstcomm);
 847
 848         /* Langevin dynamics */
 849         PR("bd-fric", ir->bd_fric);
 850         PSTEP("ld-seed", ir->ld_seed);
 851
 852         /* Energy minimization */
 853         PR("emtol", ir->em_tol);
 854         PR("emstep", ir->em_stepsize);
 855         PI("niter", ir->niter);
 856         PR("fcstep", ir->fc_stepsize);
 857         PI("nstcgsteep", ir->nstcgsteep);
 858         PI("nbfgscorr", ir->nbfgscorr);
 859
 860         /* Test particle insertion */
 861         PR("rtpi", ir->rtpi);
 862
 863         /* Output control */
 864         PI("nstxout", ir->nstxout);
 865         PI("nstvout", ir->nstvout);
 866         PI("nstfout", ir->nstfout);
 867         PI("nstlog", ir->nstlog);
 868         PI("nstcalcenergy", ir->nstcalcenergy);
 869         PI("nstenergy", ir->nstenergy);
 870         PI("nstxout-compressed", ir->nstxout_compressed);
 871         PR("compressed-x-precision", ir->x_compression_precision);
 872
 873         /* Neighborsearching parameters */
 874         PS("cutoff-scheme", ECUTSCHEME(ir->cutoff_scheme));
 875         PI("nstlist", ir->nstlist);
 876         PS("pbc", c_pbcTypeNames[ir->pbcType].c_str());
 877         PS("periodic-molecules", EBOOL(ir->bPeriodicMols));
 878         PR("verlet-buffer-tolerance", ir->verletbuf_tol);
 879         PR("rlist", ir->rlist);
 880
 881         /* Options for electrostatics and VdW */
 882         PS("coulombtype", EELTYPE(ir->coulombtype));
 883         PS("coulomb-modifier", INTMODIFIER(ir->coulomb_modifier));
 884         PR("rcoulomb-switch", ir->rcoulomb_switch);
 885         PR("rcoulomb", ir->rcoulomb);
 886         if (ir->epsilon_r != 0)
 887         {
 888             PR("epsilon-r", ir->epsilon_r);
 889         }
 890         else
 891         {
 892             PS("epsilon-r", infbuf);
 893         }
 894         if (ir->epsilon_rf != 0)
 895         {
 896             PR("epsilon-rf", ir->epsilon_rf);
 897         }
 898         else
 899         {
 900             PS("epsilon-rf", infbuf);
 901         }
 902         PS("vdw-type", EVDWTYPE(ir->vdwtype));
 903         PS("vdw-modifier", INTMODIFIER(ir->vdw_modifier));
 904         PR("rvdw-switch", ir->rvdw_switch);
 905         PR("rvdw", ir->rvdw);
 906         PS("DispCorr", EDISPCORR(ir->eDispCorr));
 907         PR("table-extension", ir->tabext);
 908
 909         PR("fourierspacing", ir->fourier_spacing);
 910         PI("fourier-nx", ir->nkx);
 911         PI("fourier-ny", ir->nky);
 912         PI("fourier-nz", ir->nkz);
 913         PI("pme-order", ir->pme_order);
 914         PR("ewald-rtol", ir->ewald_rtol);
 915         PR("ewald-rtol-lj", ir->ewald_rtol_lj);
 916         PS("lj-pme-comb-rule", ELJPMECOMBNAMES(ir->ljpme_combination_rule));
 917         PR("ewald-geometry", ir->ewald_geometry);
 918         PR("epsilon-surface", ir->epsilon_surface);
 919
 920         /* Options for weak coupling algorithms */
 921         PS("tcoupl", ETCOUPLTYPE(ir->etc));
 922         PI("nsttcouple", ir->nsttcouple);
 923         PI("nh-chain-length", ir->opts.nhchainlength);
 924         PS("print-nose-hoover-chain-variables", EBOOL(ir->bPrintNHChains));
 925
 926         PS("pcoupl", EPCOUPLTYPE(ir->epc));
 927         PS("pcoupltype", EPCOUPLTYPETYPE(ir->epct));
 928         PI("nstpcouple", ir->nstpcouple);
 929         PR("tau-p", ir->tau_p);
 930         pr_matrix(fp, indent, "compressibility", ir->compress, bMDPformat);
 931         pr_matrix(fp, indent, "ref-p", ir->ref_p, bMDPformat);
 932         PS("refcoord-scaling", EREFSCALINGTYPE(ir->refcoord_scaling));
 933
 934         if (bMDPformat)
 935         {
 936             fprintf(fp, "posres-com  = %g %g %g\n", ir->posres_com[XX], ir->posres_com[YY], ir->posres_com[ZZ]);
 937             fprintf(fp, "posres-comB = %g %g %g\n", ir->posres_comB[XX], ir->posres_comB[YY], ir->posres_comB[ZZ]);
 938         }
 939         else
 940         {
 941             pr_rvec(fp, indent, "posres-com", ir->posres_com, DIM, TRUE);
 942             pr_rvec(fp, indent, "posres-comB", ir->posres_comB, DIM, TRUE);
 943         }
 944
 945         /* QMMM */
 946         PS("QMMM", EBOOL(ir->bQMMM));
 947         fprintf(fp, "%s:\n", "qm-opts");
 948         pr_int(fp, indent, "ngQM", ir->opts.ngQM);
 949
 950         /* CONSTRAINT OPTIONS */
 951         PS("constraint-algorithm", ECONSTRTYPE(ir->eConstrAlg));
 952         PS("continuation", EBOOL(ir->bContinuation));
 953
 954         PS("Shake-SOR", EBOOL(ir->bShakeSOR));
 955         PR("shake-tol", ir->shake_tol);
 956         PI("lincs-order", ir->nProjOrder);
 957         PI("lincs-iter", ir->nLincsIter);
 958         PR("lincs-warnangle", ir->LincsWarnAngle);
 959
 960         /* Walls */
 961         PI("nwall", ir->nwall);
 962         PS("wall-type", EWALLTYPE(ir->wall_type));
 963         PR("wall-r-linpot", ir->wall_r_linpot);
 964         /* wall-atomtype */
 965         PI("wall-atomtype[0]", ir->wall_atomtype[0]);
 966         PI("wall-atomtype[1]", ir->wall_atomtype[1]);
 967         /* wall-density */
 968         PR("wall-density[0]", ir->wall_density[0]);
 969         PR("wall-density[1]", ir->wall_density[1]);
 970         PR("wall-ewald-zfac", ir->wall_ewald_zfac);
 971
 972         /* COM PULLING */
 973         PS("pull", EBOOL(ir->bPull));
 974         if (ir->bPull)
 975         {
 976             pr_pull(fp, indent, *ir->pull);
 977         }
 978
 979         /* AWH BIASING */
 980         PS("awh", EBOOL(ir->bDoAwh));
 981         if (ir->bDoAwh)
 982         {
 983             pr_awh(fp, indent, ir->awhParams);
 984         }
 985
 986         /* ENFORCED ROTATION */
 987         PS("rotation", EBOOL(ir->bRot));
 988         if (ir->bRot)
 989         {
 990             pr_rot(fp, indent, ir->rot);
 991         }
 992
 993         /* INTERACTIVE MD */
 994         PS("interactiveMD", EBOOL(ir->bIMD));
 995         if (ir->bIMD)
 996         {
 997             pr_imd(fp, indent, ir->imd);
 998         }
 999
1000         /* NMR refinement stuff */
1001         PS("disre", EDISRETYPE(ir->eDisre));
1002         PS("disre-weighting", EDISREWEIGHTING(ir->eDisreWeighting));
1003         PS("disre-mixed", EBOOL(ir->bDisreMixed));
1004         PR("dr-fc", ir->dr_fc);
1005         PR("dr-tau", ir->dr_tau);
1006         PR("nstdisreout", ir->nstdisreout);
1007
1008         PR("orire-fc", ir->orires_fc);
1009         PR("orire-tau", ir->orires_tau);
1010         PR("nstorireout", ir->nstorireout);
1011
1012         /* FREE ENERGY VARIABLES */
1013         PS("free-energy", EFEPTYPE(ir->efep));
1014         if (ir->efep != efepNO || ir->bSimTemp)
1015         {
1016             pr_fepvals(fp, indent, ir->fepvals, bMDPformat);
1017         }
1018         if (ir->bExpanded)
1019         {
1020             pr_expandedvals(fp, indent, ir->expandedvals, ir->fepvals->n_lambda);
1021         }
1022
1023         /* NON-equilibrium MD stuff */
1024         PR("cos-acceleration", ir->cos_accel);
1025         pr_matrix(fp, indent, "deform", ir->deform, bMDPformat);
1026
1027         /* SIMULATED TEMPERING */
1028         PS("simulated-tempering", EBOOL(ir->bSimTemp));
1029         if (ir->bSimTemp)
1030         {
1031             pr_simtempvals(fp, indent, ir->simtempvals, ir->fepvals->n_lambda);
1032         }
1033
1034         /* ION/WATER SWAPPING FOR COMPUTATIONAL ELECTROPHYSIOLOGY */
1035         PS("swapcoords", ESWAPTYPE(ir->eSwapCoords));
1036         if (ir->eSwapCoords != eswapNO)
1037         {
1038             pr_swap(fp, indent, ir->swap);
1039         }
1040
1041         /* USER-DEFINED THINGIES */
1042         PI("userint1", ir->userint1);
1043         PI("userint2", ir->userint2);
1044         PI("userint3", ir->userint3);
1045         PI("userint4", ir->userint4);
1046         PR("userreal1", ir->userreal1);
1047         PR("userreal2", ir->userreal2);
1048         PR("userreal3", ir->userreal3);
1049         PR("userreal4", ir->userreal4);
1050
1051         if (!bMDPformat)
1052         {
1053             gmx::TextWriter writer(fp);
1054             writer.wrapperSettings().setIndent(indent);
1055             gmx::dumpKeyValueTree(&writer, *ir->params);
1056         }
1057
1058         pr_grp_opts(fp, indent, "grpopts", &(ir->opts), bMDPformat);
1059     }
1060 }
1061 #undef PS
1062 #undef PR
1063 #undef PI
1064
1065 static void cmp_grpopts(FILE* fp, const t_grpopts* opt1, const t_grpopts* opt2, real ftol, real abstol)
1066 {
1067     int  i, j;
1068     char buf1[256], buf2[256];
1069
1070     cmp_int(fp, "inputrec->grpopts.ngtc", -1, opt1->ngtc, opt2->ngtc);
1071     cmp_int(fp, "inputrec->grpopts.ngacc", -1, opt1->ngacc, opt2->ngacc);
1072     cmp_int(fp, "inputrec->grpopts.ngfrz", -1, opt1->ngfrz, opt2->ngfrz);
1073     cmp_int(fp, "inputrec->grpopts.ngener", -1, opt1->ngener, opt2->ngener);
1074     for (i = 0; (i < std::min(opt1->ngtc, opt2->ngtc)); i++)
1075     {
1076         cmp_real(fp, "inputrec->grpopts.nrdf", i, opt1->nrdf[i], opt2->nrdf[i], ftol, abstol);
1077         cmp_real(fp, "inputrec->grpopts.ref_t", i, opt1->ref_t[i], opt2->ref_t[i], ftol, abstol);
1078         cmp_real(fp, "inputrec->grpopts.tau_t", i, opt1->tau_t[i], opt2->tau_t[i], ftol, abstol);
1079         cmp_int(fp, "inputrec->grpopts.annealing", i, opt1->annealing[i], opt2->annealing[i]);
1080         cmp_int(fp, "inputrec->grpopts.anneal_npoints", i, opt1->anneal_npoints[i], opt2->anneal_npoints[i]);
1081         if (opt1->anneal_npoints[i] == opt2->anneal_npoints[i])
1082         {
1083             sprintf(buf1, "inputrec->grpopts.anneal_time[%d]", i);
1084             sprintf(buf2, "inputrec->grpopts.anneal_temp[%d]", i);
1085             for (j = 0; j < opt1->anneal_npoints[i]; j++)
1086             {
1087                 cmp_real(fp, buf1, j, opt1->anneal_time[i][j], opt2->anneal_time[i][j], ftol, abstol);
1088                 cmp_real(fp, buf2, j, opt1->anneal_temp[i][j], opt2->anneal_temp[i][j], ftol, abstol);
1089             }
1090         }
1091     }
1092     if (opt1->ngener == opt2->ngener)
1093     {
1094         for (i = 0; i < opt1->ngener; i++)
1095         {
1096             for (j = i; j < opt1->ngener; j++)
1097             {
1098                 sprintf(buf1, "inputrec->grpopts.egp_flags[%d]", i);
1099                 cmp_int(fp, buf1, j, opt1->egp_flags[opt1->ngener * i + j], opt2->egp_flags[opt1->ngener * i + j]);
1100             }
1101         }
1102     }
1103     for (i = 0; (i < std::min(opt1->ngacc, opt2->ngacc)); i++)
1104     {
1105         cmp_rvec(fp, "inputrec->grpopts.acc", i, opt1->acc[i], opt2->acc[i], ftol, abstol);
1106     }
1107     for (i = 0; (i < std::min(opt1->ngfrz, opt2->ngfrz)); i++)
1108     {
1109         cmp_ivec(fp, "inputrec->grpopts.nFreeze", i, opt1->nFreeze[i], opt2->nFreeze[i]);
1110     }
1111 }
1112
1113 static void cmp_pull(FILE* fp)
1114 {
1115     fprintf(fp,
1116             "WARNING: Both files use COM pulling, but comparing of the pull struct is not "
1117             "implemented (yet). The pull parameters could be the same or different.\n");
1118 }
1119
1120 static void cmp_awhDimParams(FILE*                    fp,
1121                              const gmx::AwhDimParams* dimp1,
1122                              const gmx::AwhDimParams* dimp2,
1123                              int                      dimIndex,
1124                              real                     ftol,
1125                              real                     abstol)
1126 {
1127     /* Note that we have double index here, but the compare functions only
1128      * support one index, so here we only print the dim index and not the bias.
1129      */
1130     cmp_int(fp, "inputrec.awhParams->bias?->dim->coord_index", dimIndex, dimp1->coordIndex, dimp2->coordIndex);
1131     cmp_double(fp, "inputrec->awhParams->bias?->dim->period", dimIndex, dimp1->period, dimp2->period, ftol, abstol);
1132     cmp_double(fp, "inputrec->awhParams->bias?->dim->diffusion", dimIndex, dimp1->diffusion, dimp2->diffusion, ftol, abstol);
1133     cmp_double(fp, "inputrec->awhParams->bias?->dim->origin", dimIndex, dimp1->origin, dimp2->origin, ftol, abstol);
1134     cmp_double(fp, "inputrec->awhParams->bias?->dim->end", dimIndex, dimp1->end, dimp2->end, ftol, abstol);
1135     cmp_double(fp,
1136                "inputrec->awhParams->bias?->dim->coord_value_init",
1137                dimIndex,
1138                dimp1->coordValueInit,
1139                dimp2->coordValueInit,
1140                ftol,
1141                abstol);
1142     cmp_double(fp,
1143                "inputrec->awhParams->bias?->dim->coverDiameter",
1144                dimIndex,
1145                dimp1->coverDiameter,
1146                dimp2->coverDiameter,
1147                ftol,
1148                abstol);
1149 }
1150
1151 static void cmp_awhBiasParams(FILE*                     fp,
1152                               const gmx::AwhBiasParams* bias1,
1153                               const gmx::AwhBiasParams* bias2,
1154                               int                       biasIndex,
1155                               real                      ftol,
1156                               real                      abstol)
1157 {
1158     cmp_int(fp, "inputrec->awhParams->ndim", biasIndex, bias1->ndim, bias2->ndim);
1159     cmp_int(fp, "inputrec->awhParams->biaseTarget", biasIndex, bias1->eTarget, bias2->eTarget);
1160     cmp_double(fp,
1161                "inputrec->awhParams->biastargetBetaScaling",
1162                biasIndex,
1163                bias1->targetBetaScaling,
1164                bias2->targetBetaScaling,
1165                ftol,
1166                abstol);
1167     cmp_double(fp,
1168                "inputrec->awhParams->biastargetCutoff",
1169                biasIndex,
1170                bias1->targetCutoff,
1171                bias2->targetCutoff,
1172                ftol,
1173                abstol);
1174     cmp_int(fp, "inputrec->awhParams->biaseGrowth", biasIndex, bias1->eGrowth, bias2->eGrowth);
1175     cmp_bool(fp, "inputrec->awhParams->biasbUserData", biasIndex, bias1->bUserData != 0, bias2->bUserData != 0);
1176     cmp_double(fp,
1177                "inputrec->awhParams->biaserror_initial",
1178                biasIndex,
1179                bias1->errorInitial,
1180                bias2->errorInitial,
1181                ftol,
1182                abstol);
1183     cmp_int(fp, "inputrec->awhParams->biasShareGroup", biasIndex, bias1->shareGroup, bias2->shareGroup);
1184
1185     for (int dim = 0; dim < std::min(bias1->ndim, bias2->ndim); dim++)
1186     {
1187         cmp_awhDimParams(fp, &bias1->dimParams[dim], &bias2->dimParams[dim], dim, ftol, abstol);
1188     }
1189 }
1190
1191 static void cmp_awhParams(FILE* fp, const gmx::AwhParams* awh1, const gmx::AwhParams* awh2, real ftol, real abstol)
1192 {
1193     cmp_int(fp, "inputrec->awhParams->nbias", -1, awh1->numBias, awh2->numBias);
1194     cmp_int64(fp, "inputrec->awhParams->seed", awh1->seed, awh2->seed);
1195     cmp_int(fp, "inputrec->awhParams->nstout", -1, awh1->nstOut, awh2->nstOut);
1196     cmp_int(fp, "inputrec->awhParams->nstsample_coord", -1, awh1->nstSampleCoord, awh2->nstSampleCoord);
1197     cmp_int(fp,
1198             "inputrec->awhParams->nsamples_update_free_energy",
1199             -1,
1200             awh1->numSamplesUpdateFreeEnergy,
1201             awh2->numSamplesUpdateFreeEnergy);
1202     cmp_int(fp, "inputrec->awhParams->ePotential", -1, awh1->ePotential, awh2->ePotential);
1203     cmp_bool(fp, "inputrec->awhParams->shareBiasMultisim", -1, awh1->shareBiasMultisim, awh2->shareBiasMultisim);
1204
1205     if (awh1->numBias == awh2->numBias)
1206     {
1207         for (int bias = 0; bias < awh1->numBias; bias++)
1208         {
1209             cmp_awhBiasParams(fp, &awh1->awhBiasParams[bias], &awh2->awhBiasParams[bias], bias, ftol, abstol);
1210         }
1211     }
1212 }
1213
1214 static void cmp_simtempvals(FILE*            fp,
1215                             const t_simtemp* simtemp1,
1216                             const t_simtemp* simtemp2,
1217                             int              n_lambda,
1218                             real             ftol,
1219                             real             abstol)
1220 {
1221     int i;
1222     cmp_int(fp, "inputrec->simtempvals->eSimTempScale", -1, simtemp1->eSimTempScale, simtemp2->eSimTempScale);
1223     cmp_real(fp, "inputrec->simtempvals->simtemp_high", -1, simtemp1->simtemp_high, simtemp2->simtemp_high, ftol, abstol);
1224     cmp_real(fp, "inputrec->simtempvals->simtemp_low", -1, simtemp1->simtemp_low, simtemp2->simtemp_low, ftol, abstol);
1225     for (i = 0; i < n_lambda; i++)
1226     {
1227         cmp_real(fp,
1228                  "inputrec->simtempvals->temperatures",
1229                  -1,
1230                  simtemp1->temperatures[i],
1231                  simtemp2->temperatures[i],
1232                  ftol,
1233                  abstol);
1234     }
1235 }
1236
1237 static void cmp_expandedvals(FILE*             fp,
1238                              const t_expanded* expand1,
1239                              const t_expanded* expand2,
1240                              int               n_lambda,
1241                              real              ftol,
1242                              real              abstol)
1243 {
1244     int i;
1245
1246     cmp_bool(fp, "inputrec->fepvals->bInit_weights", -1, expand1->bInit_weights, expand2->bInit_weights);
1247     cmp_bool(fp, "inputrec->fepvals->bWLoneovert", -1, expand1->bWLoneovert, expand2->bWLoneovert);
1248
1249     for (i = 0; i < n_lambda; i++)
1250     {
1251         cmp_real(fp,
1252                  "inputrec->expandedvals->init_lambda_weights",
1253                  -1,
1254                  expand1->init_lambda_weights[i],
1255                  expand2->init_lambda_weights[i],
1256                  ftol,
1257                  abstol);
1258     }
1259
1260     cmp_int(fp, "inputrec->expandedvals->lambda-stats", -1, expand1->elamstats, expand2->elamstats);
1261     cmp_int(fp, "inputrec->expandedvals->lambda-mc-move", -1, expand1->elmcmove, expand2->elmcmove);
1262     cmp_int(fp, "inputrec->expandedvals->lmc-repeats", -1, expand1->lmc_repeats, expand2->lmc_repeats);
1263     cmp_int(fp, "inputrec->expandedvals->lmc-gibbsdelta", -1, expand1->gibbsdeltalam, expand2->gibbsdeltalam);
1264     cmp_int(fp, "inputrec->expandedvals->lmc-forced-nstart", -1, expand1->lmc_forced_nstart, expand2->lmc_forced_nstart);
1265     cmp_int(fp, "inputrec->expandedvals->lambda-weights-equil", -1, expand1->elmceq, expand2->elmceq);
1266     cmp_int(fp,
1267             "inputrec->expandedvals->,weight-equil-number-all-lambda",
1268             -1,
1269             expand1->equil_n_at_lam,
1270             expand2->equil_n_at_lam);
1271     cmp_int(fp, "inputrec->expandedvals->weight-equil-number-samples", -1, expand1->equil_samples, expand2->equil_samples);
1272     cmp_int(fp, "inputrec->expandedvals->weight-equil-number-steps", -1, expand1->equil_steps, expand2->equil_steps);
1273     cmp_real(fp,
1274              "inputrec->expandedvals->weight-equil-wl-delta",
1275              -1,
1276              expand1->equil_wl_delta,
1277              expand2->equil_wl_delta,
1278              ftol,
1279              abstol);
1280     cmp_real(fp,
1281              "inputrec->expandedvals->weight-equil-count-ratio",
1282              -1,
1283              expand1->equil_ratio,
1284              expand2->equil_ratio,
1285              ftol,
1286              abstol);
1287     cmp_bool(fp,
1288              "inputrec->expandedvals->symmetrized-transition-matrix",
1289              -1,
1290              expand1->bSymmetrizedTMatrix,
1291              expand2->bSymmetrizedTMatrix);
1292     cmp_int(fp, "inputrec->expandedvals->nstTij", -1, expand1->nstTij, expand2->nstTij);
1293     cmp_int(fp, "inputrec->expandedvals->mininum-var-min", -1, expand1->minvarmin, expand2->minvarmin); /*default is reasonable */
1294     cmp_int(fp, "inputrec->expandedvals->weight-c-range", -1, expand1->c_range, expand2->c_range); /* default is just C=0 */
1295     cmp_real(fp, "inputrec->expandedvals->wl-scale", -1, expand1->wl_scale, expand2->wl_scale, ftol, abstol);
1296     cmp_real(fp, "inputrec->expandedvals->init-wl-delta", -1, expand1->init_wl_delta, expand2->init_wl_delta, ftol, abstol);
1297     cmp_real(fp, "inputrec->expandedvals->wl-ratio", -1, expand1->wl_ratio, expand2->wl_ratio, ftol, abstol);
1298     cmp_int(fp, "inputrec->expandedvals->nstexpanded", -1, expand1->nstexpanded, expand2->nstexpanded);
1299     cmp_int(fp, "inputrec->expandedvals->lmc-seed", -1, expand1->lmc_seed, expand2->lmc_seed);
1300     cmp_real(fp, "inputrec->expandedvals->mc-temperature", -1, expand1->mc_temp, expand2->mc_temp, ftol, abstol);
1301 }
1302
1303 static void cmp_fepvals(FILE* fp, const t_lambda* fep1, const t_lambda* fep2, real ftol, real abstol)
1304 {
1305     int i, j;
1306     cmp_int(fp, "inputrec->nstdhdl", -1, fep1->nstdhdl, fep2->nstdhdl);
1307     cmp_double(fp, "inputrec->fepvals->init_fep_state", -1, fep1->init_fep_state, fep2->init_fep_state, ftol, abstol);
1308     cmp_double(fp, "inputrec->fepvals->delta_lambda", -1, fep1->delta_lambda, fep2->delta_lambda, ftol, abstol);
1309     cmp_int(fp, "inputrec->fepvals->n_lambda", -1, fep1->n_lambda, fep2->n_lambda);
1310     for (i = 0; i < efptNR; i++)
1311     {
1312         for (j = 0; j < std::min(fep1->n_lambda, fep2->n_lambda); j++)
1313         {
1314             cmp_double(fp,
1315                        "inputrec->fepvals->all_lambda",
1316                        -1,
1317                        fep1->all_lambda[i][j],
1318                        fep2->all_lambda[i][j],
1319                        ftol,
1320                        abstol);
1321         }
1322     }
1323     cmp_int(fp, "inputrec->fepvals->lambda_neighbors", 1, fep1->lambda_neighbors, fep2->lambda_neighbors);
1324     cmp_real(fp, "inputrec->fepvals->sc_alpha", -1, fep1->sc_alpha, fep2->sc_alpha, ftol, abstol);
1325     cmp_int(fp, "inputrec->fepvals->sc_power", -1, fep1->sc_power, fep2->sc_power);
1326     cmp_real(fp, "inputrec->fepvals->sc_r_power", -1, fep1->sc_r_power, fep2->sc_r_power, ftol, abstol);
1327     cmp_real(fp, "inputrec->fepvals->sc_sigma", -1, fep1->sc_sigma, fep2->sc_sigma, ftol, abstol);
1328     cmp_int(fp, "inputrec->fepvals->edHdLPrintEnergy", -1, fep1->edHdLPrintEnergy, fep1->edHdLPrintEnergy);
1329     cmp_bool(fp, "inputrec->fepvals->bScCoul", -1, fep1->bScCoul, fep1->bScCoul);
1330     cmp_int(fp, "inputrec->separate_dhdl_file", -1, fep1->separate_dhdl_file, fep2->separate_dhdl_file);
1331     cmp_int(fp, "inputrec->dhdl_derivatives", -1, fep1->dhdl_derivatives, fep2->dhdl_derivatives);
1332     cmp_int(fp, "inputrec->dh_hist_size", -1, fep1->dh_hist_size, fep2->dh_hist_size);
1333     cmp_double(fp, "inputrec->dh_hist_spacing", -1, fep1->dh_hist_spacing, fep2->dh_hist_spacing, ftol, abstol);
1334 }
1335
1336 void cmp_inputrec(FILE* fp, const t_inputrec* ir1, const t_inputrec* ir2, real ftol, real abstol)
1337 {
1338     fprintf(fp, "comparing inputrec\n");
1339
1340     /* gcc 2.96 doesnt like these defines at all, but issues a huge list
1341      * of warnings. Maybe it will change in future versions, but for the
1342      * moment I've spelled them out instead. /EL 000820
1343      * #define CIB(s) cmp_int(fp,"inputrec->"#s,0,ir1->##s,ir2->##s)
1344      * #define CII(s) cmp_int(fp,"inputrec->"#s,0,ir1->##s,ir2->##s)
1345      * #define CIR(s) cmp_real(fp,"inputrec->"#s,0,ir1->##s,ir2->##s,ftol)
1346      */
1347     cmp_int(fp, "inputrec->eI", -1, ir1->eI, ir2->eI);
1348     cmp_int64(fp, "inputrec->nsteps", ir1->nsteps, ir2->nsteps);
1349     cmp_int64(fp, "inputrec->init_step", ir1->init_step, ir2->init_step);
1350     cmp_int(fp, "inputrec->simulation_part", -1, ir1->simulation_part, ir2->simulation_part);
1351     cmp_int(fp, "inputrec->mts", -1, static_cast<int>(ir1->useMts), static_cast<int>(ir2->useMts));
1352     if (ir1->useMts && ir2->useMts)
1353     {
1354         cmp_int(fp, "inputrec->mts-levels", -1, ir1->mtsLevels.size(), ir2->mtsLevels.size());
1355         cmp_int(fp,
1356                 "inputrec->mts-level2-forces",
1357                 -1,
1358                 ir1->mtsLevels[1].forceGroups.to_ulong(),
1359                 ir2->mtsLevels[1].forceGroups.to_ulong());
1360         cmp_int(fp,
1361                 "inputrec->mts-level2-factor",
1362                 -1,
1363                 ir1->mtsLevels[1].stepFactor,
1364                 ir2->mtsLevels[1].stepFactor);
1365     }
1366     cmp_int(fp, "inputrec->pbcType", -1, static_cast<int>(ir1->pbcType), static_cast<int>(ir2->pbcType));
1367     cmp_bool(fp, "inputrec->bPeriodicMols", -1, ir1->bPeriodicMols, ir2->bPeriodicMols);
1368     cmp_int(fp, "inputrec->cutoff_scheme", -1, ir1->cutoff_scheme, ir2->cutoff_scheme);
1369     cmp_int(fp, "inputrec->nstlist", -1, ir1->nstlist, ir2->nstlist);
1370     cmp_int(fp, "inputrec->nstcomm", -1, ir1->nstcomm, ir2->nstcomm);
1371     cmp_int(fp, "inputrec->comm_mode", -1, ir1->comm_mode, ir2->comm_mode);
1372     cmp_int(fp, "inputrec->nstlog", -1, ir1->nstlog, ir2->nstlog);
1373     cmp_int(fp, "inputrec->nstxout", -1, ir1->nstxout, ir2->nstxout);
1374     cmp_int(fp, "inputrec->nstvout", -1, ir1->nstvout, ir2->nstvout);
1375     cmp_int(fp, "inputrec->nstfout", -1, ir1->nstfout, ir2->nstfout);
1376     cmp_int(fp, "inputrec->nstcalcenergy", -1, ir1->nstcalcenergy, ir2->nstcalcenergy);
1377     cmp_int(fp, "inputrec->nstenergy", -1, ir1->nstenergy, ir2->nstenergy);
1378     cmp_int(fp, "inputrec->nstxout_compressed", -1, ir1->nstxout_compressed, ir2->nstxout_compressed);
1379     cmp_double(fp, "inputrec->init_t", -1, ir1->init_t, ir2->init_t, ftol, abstol);
1380     cmp_double(fp, "inputrec->delta_t", -1, ir1->delta_t, ir2->delta_t, ftol, abstol);
1381     cmp_real(fp,
1382              "inputrec->x_compression_precision",
1383              -1,
1384              ir1->x_compression_precision,
1385              ir2->x_compression_precision,
1386              ftol,
1387              abstol);
1388     cmp_real(fp, "inputrec->fourierspacing", -1, ir1->fourier_spacing, ir2->fourier_spacing, ftol, abstol);
1389     cmp_int(fp, "inputrec->nkx", -1, ir1->nkx, ir2->nkx);
1390     cmp_int(fp, "inputrec->nky", -1, ir1->nky, ir2->nky);
1391     cmp_int(fp, "inputrec->nkz", -1, ir1->nkz, ir2->nkz);
1392     cmp_int(fp, "inputrec->pme_order", -1, ir1->pme_order, ir2->pme_order);
1393     cmp_real(fp, "inputrec->ewald_rtol", -1, ir1->ewald_rtol, ir2->ewald_rtol, ftol, abstol);
1394     cmp_int(fp, "inputrec->ewald_geometry", -1, ir1->ewald_geometry, ir2->ewald_geometry);
1395     cmp_real(fp, "inputrec->epsilon_surface", -1, ir1->epsilon_surface, ir2->epsilon_surface, ftol, abstol);
1396     cmp_int(fp,
1397             "inputrec->bContinuation",
1398             -1,
1399             static_cast<int>(ir1->bContinuation),
1400             static_cast<int>(ir2->bContinuation));
1401     cmp_int(fp, "inputrec->bShakeSOR", -1, static_cast<int>(ir1->bShakeSOR), static_cast<int>(ir2->bShakeSOR));
1402     cmp_int(fp, "inputrec->etc", -1, ir1->etc, ir2->etc);
1403     cmp_int(fp,
1404             "inputrec->bPrintNHChains",
1405             -1,
1406             static_cast<int>(ir1->bPrintNHChains),
1407             static_cast<int>(ir2->bPrintNHChains));
1408     cmp_int(fp, "inputrec->epc", -1, ir1->epc, ir2->epc);
1409     cmp_int(fp, "inputrec->epct", -1, ir1->epct, ir2->epct);
1410     cmp_real(fp, "inputrec->tau_p", -1, ir1->tau_p, ir2->tau_p, ftol, abstol);
1411     cmp_rvec(fp, "inputrec->ref_p(x)", -1, ir1->ref_p[XX], ir2->ref_p[XX], ftol, abstol);
1412     cmp_rvec(fp, "inputrec->ref_p(y)", -1, ir1->ref_p[YY], ir2->ref_p[YY], ftol, abstol);
1413     cmp_rvec(fp, "inputrec->ref_p(z)", -1, ir1->ref_p[ZZ], ir2->ref_p[ZZ], ftol, abstol);
1414     cmp_rvec(fp, "inputrec->compress(x)", -1, ir1->compress[XX], ir2->compress[XX], ftol, abstol);
1415     cmp_rvec(fp, "inputrec->compress(y)", -1, ir1->compress[YY], ir2->compress[YY], ftol, abstol);
1416     cmp_rvec(fp, "inputrec->compress(z)", -1, ir1->compress[ZZ], ir2->compress[ZZ], ftol, abstol);
1417     cmp_int(fp, "refcoord_scaling", -1, ir1->refcoord_scaling, ir2->refcoord_scaling);
1418     cmp_rvec(fp, "inputrec->posres_com", -1, ir1->posres_com, ir2->posres_com, ftol, abstol);
1419     cmp_rvec(fp, "inputrec->posres_comB", -1, ir1->posres_comB, ir2->posres_comB, ftol, abstol);
1420     cmp_real(fp, "inputrec->verletbuf_tol", -1, ir1->verletbuf_tol, ir2->verletbuf_tol, ftol, abstol);
1421     cmp_real(fp, "inputrec->rlist", -1, ir1->rlist, ir2->rlist, ftol, abstol);
1422     cmp_real(fp, "inputrec->rtpi", -1, ir1->rtpi, ir2->rtpi, ftol, abstol);
1423     cmp_int(fp, "inputrec->coulombtype", -1, ir1->coulombtype, ir2->coulombtype);
1424     cmp_int(fp, "inputrec->coulomb_modifier", -1, ir1->coulomb_modifier, ir2->coulomb_modifier);
1425     cmp_real(fp, "inputrec->rcoulomb_switch", -1, ir1->rcoulomb_switch, ir2->rcoulomb_switch, ftol, abstol);
1426     cmp_real(fp, "inputrec->rcoulomb", -1, ir1->rcoulomb, ir2->rcoulomb, ftol, abstol);
1427     cmp_int(fp, "inputrec->vdwtype", -1, ir1->vdwtype, ir2->vdwtype);
1428     cmp_int(fp, "inputrec->vdw_modifier", -1, ir1->vdw_modifier, ir2->vdw_modifier);
1429     cmp_real(fp, "inputrec->rvdw_switch", -1, ir1->rvdw_switch, ir2->rvdw_switch, ftol, abstol);
1430     cmp_real(fp, "inputrec->rvdw", -1, ir1->rvdw, ir2->rvdw, ftol, abstol);
1431     cmp_real(fp, "inputrec->epsilon_r", -1, ir1->epsilon_r, ir2->epsilon_r, ftol, abstol);
1432     cmp_real(fp, "inputrec->epsilon_rf", -1, ir1->epsilon_rf, ir2->epsilon_rf, ftol, abstol);
1433     cmp_real(fp, "inputrec->tabext", -1, ir1->tabext, ir2->tabext, ftol, abstol);
1434
1435     cmp_int(fp, "inputrec->eDispCorr", -1, ir1->eDispCorr, ir2->eDispCorr);
1436     cmp_real(fp, "inputrec->shake_tol", -1, ir1->shake_tol, ir2->shake_tol, ftol, abstol);
1437     cmp_int(fp, "inputrec->efep", -1, ir1->efep, ir2->efep);
1438     cmp_fepvals(fp, ir1->fepvals, ir2->fepvals, ftol, abstol);
1439     cmp_int(fp, "inputrec->bSimTemp", -1, static_cast<int>(ir1->bSimTemp), static_cast<int>(ir2->bSimTemp));
1440     if ((ir1->bSimTemp == ir2->bSimTemp) && (ir1->bSimTemp))
1441     {
1442         cmp_simtempvals(fp,
1443                         ir1->simtempvals,
1444                         ir2->simtempvals,
1445                         std::min(ir1->fepvals->n_lambda, ir2->fepvals->n_lambda),
1446                         ftol,
1447                         abstol);
1448     }
1449     cmp_int(fp, "inputrec->bExpanded", -1, static_cast<int>(ir1->bExpanded), static_cast<int>(ir2->bExpanded));
1450     if ((ir1->bExpanded == ir2->bExpanded) && (ir1->bExpanded))
1451     {
1452         cmp_expandedvals(fp,
1453                          ir1->expandedvals,
1454                          ir2->expandedvals,
1455                          std::min(ir1->fepvals->n_lambda, ir2->fepvals->n_lambda),
1456                          ftol,
1457                          abstol);
1458     }
1459     cmp_int(fp, "inputrec->nwall", -1, ir1->nwall, ir2->nwall);
1460     cmp_int(fp, "inputrec->wall_type", -1, ir1->wall_type, ir2->wall_type);
1461     cmp_int(fp, "inputrec->wall_atomtype[0]", -1, ir1->wall_atomtype[0], ir2->wall_atomtype[0]);
1462     cmp_int(fp, "inputrec->wall_atomtype[1]", -1, ir1->wall_atomtype[1], ir2->wall_atomtype[1]);
1463     cmp_real(fp, "inputrec->wall_density[0]", -1, ir1->wall_density[0], ir2->wall_density[0], ftol, abstol);
1464     cmp_real(fp, "inputrec->wall_density[1]", -1, ir1->wall_density[1], ir2->wall_density[1], ftol, abstol);
1465     cmp_real(fp, "inputrec->wall_ewald_zfac", -1, ir1->wall_ewald_zfac, ir2->wall_ewald_zfac, ftol, abstol);
1466
1467     cmp_bool(fp, "inputrec->bPull", -1, ir1->bPull, ir2->bPull);
1468     if (ir1->bPull && ir2->bPull)
1469     {
1470         cmp_pull(fp);
1471     }
1472
1473     cmp_bool(fp, "inputrec->bDoAwh", -1, ir1->bDoAwh, ir2->bDoAwh);
1474     if (ir1->bDoAwh && ir2->bDoAwh)
1475     {
1476         cmp_awhParams(fp, ir1->awhParams, ir2->awhParams, ftol, abstol);
1477     }
1478
1479     cmp_int(fp, "inputrec->eDisre", -1, ir1->eDisre, ir2->eDisre);
1480     cmp_real(fp, "inputrec->dr_fc", -1, ir1->dr_fc, ir2->dr_fc, ftol, abstol);
1481     cmp_int(fp, "inputrec->eDisreWeighting", -1, ir1->eDisreWeighting, ir2->eDisreWeighting);
1482     cmp_int(fp, "inputrec->bDisreMixed", -1, static_cast<int>(ir1->bDisreMixed), static_cast<int>(ir2->bDisreMixed));
1483     cmp_int(fp, "inputrec->nstdisreout", -1, ir1->nstdisreout, ir2->nstdisreout);
1484     cmp_real(fp, "inputrec->dr_tau", -1, ir1->dr_tau, ir2->dr_tau, ftol, abstol);
1485     cmp_real(fp, "inputrec->orires_fc", -1, ir1->orires_fc, ir2->orires_fc, ftol, abstol);
1486     cmp_real(fp, "inputrec->orires_tau", -1, ir1->orires_tau, ir2->orires_tau, ftol, abstol);
1487     cmp_int(fp, "inputrec->nstorireout", -1, ir1->nstorireout, ir2->nstorireout);
1488     cmp_real(fp, "inputrec->em_stepsize", -1, ir1->em_stepsize, ir2->em_stepsize, ftol, abstol);
1489     cmp_real(fp, "inputrec->em_tol", -1, ir1->em_tol, ir2->em_tol, ftol, abstol);
1490     cmp_int(fp, "inputrec->niter", -1, ir1->niter, ir2->niter);
1491     cmp_real(fp, "inputrec->fc_stepsize", -1, ir1->fc_stepsize, ir2->fc_stepsize, ftol, abstol);
1492     cmp_int(fp, "inputrec->nstcgsteep", -1, ir1->nstcgsteep, ir2->nstcgsteep);
1493     cmp_int(fp, "inputrec->nbfgscorr", 0, ir1->nbfgscorr, ir2->nbfgscorr);
1494     cmp_int(fp, "inputrec->eConstrAlg", -1, ir1->eConstrAlg, ir2->eConstrAlg);
1495     cmp_int(fp, "inputrec->nProjOrder", -1, ir1->nProjOrder, ir2->nProjOrder);
1496     cmp_real(fp, "inputrec->LincsWarnAngle", -1, ir1->LincsWarnAngle, ir2->LincsWarnAngle, ftol, abstol);
1497     cmp_int(fp, "inputrec->nLincsIter", -1, ir1->nLincsIter, ir2->nLincsIter);
1498     cmp_real(fp, "inputrec->bd_fric", -1, ir1->bd_fric, ir2->bd_fric, ftol, abstol);
1499     cmp_int64(fp, "inputrec->ld_seed", ir1->ld_seed, ir2->ld_seed);
1500     cmp_real(fp, "inputrec->cos_accel", -1, ir1->cos_accel, ir2->cos_accel, ftol, abstol);
1501     cmp_rvec(fp, "inputrec->deform(a)", -1, ir1->deform[XX], ir2->deform[XX], ftol, abstol);
1502     cmp_rvec(fp, "inputrec->deform(b)", -1, ir1->deform[YY], ir2->deform[YY], ftol, abstol);
1503     cmp_rvec(fp, "inputrec->deform(c)", -1, ir1->deform[ZZ], ir2->deform[ZZ], ftol, abstol);
1504
1505
1506     cmp_int(fp, "inputrec->userint1", -1, ir1->userint1, ir2->userint1);
1507     cmp_int(fp, "inputrec->userint2", -1, ir1->userint2, ir2->userint2);
1508     cmp_int(fp, "inputrec->userint3", -1, ir1->userint3, ir2->userint3);
1509     cmp_int(fp, "inputrec->userint4", -1, ir1->userint4, ir2->userint4);
1510     cmp_real(fp, "inputrec->userreal1", -1, ir1->userreal1, ir2->userreal1, ftol, abstol);
1511     cmp_real(fp, "inputrec->userreal2", -1, ir1->userreal2, ir2->userreal2, ftol, abstol);
1512     cmp_real(fp, "inputrec->userreal3", -1, ir1->userreal3, ir2->userreal3, ftol, abstol);
1513     cmp_real(fp, "inputrec->userreal4", -1, ir1->userreal4, ir2->userreal4, ftol, abstol);
1514     cmp_grpopts(fp, &(ir1->opts), &(ir2->opts), ftol, abstol);
1515     gmx::TextWriter writer(fp);
1516     gmx::compareKeyValueTrees(&writer, *ir1->params, *ir2->params, ftol, abstol);
1517 }
1518
1519 void comp_pull_AB(FILE* fp, const pull_params_t& pull, real ftol, real abstol)
1520 {
1521     for (int i = 0; i < pull.ncoord; i++)
1522     {
1523         fprintf(fp, "comparing pull coord %d\n", i);
1524         cmp_real(fp, "pull-coord->k", -1, pull.coord[i].k, pull.coord[i].kB, ftol, abstol);
1525     }
1526 }
1527
1528 gmx_bool inputrecDeform(const t_inputrec* ir)
1529 {
1530     return (ir->deform[XX][XX] != 0 || ir->deform[YY][YY] != 0 || ir->deform[ZZ][ZZ] != 0
1531             || ir->deform[YY][XX] != 0 || ir->deform[ZZ][XX] != 0 || ir->deform[ZZ][YY] != 0);
1532 }
1533
1534 gmx_bool inputrecDynamicBox(const t_inputrec* ir)
1535 {
1536     return (ir->epc != epcNO || ir->eI == eiTPI || inputrecDeform(ir));
1537 }
1538
1539 gmx_bool inputrecPreserveShape(const t_inputrec* ir)
1540 {
1541     return (ir->epc != epcNO && ir->deform[XX][XX] == 0
1542             && (ir->epct == epctISOTROPIC || ir->epct == epctSEMIISOTROPIC));
1543 }
1544
1545 gmx_bool inputrecNeedMutot(const t_inputrec* ir)
1546 {
1547     return ((ir->coulombtype == eelEWALD || EEL_PME(ir->coulombtype))
1548             && (ir->ewald_geometry == eewg3DC || ir->epsilon_surface != 0));
1549 }
1550
1551 gmx_bool inputrecExclForces(const t_inputrec* ir)
1552 {
1553     return (EEL_FULL(ir->coulombtype) || (EEL_RF(ir->coulombtype)));
1554 }
1555
1556 gmx_bool inputrecNptTrotter(const t_inputrec* ir)
1557 {
1558     return (((ir->eI == eiVV) || (ir->eI == eiVVAK)) && (ir->epc == epcMTTK) && (ir->etc == etcNOSEHOOVER));
1559 }
1560
1561 gmx_bool inputrecNvtTrotter(const t_inputrec* ir)
1562 {
1563     return (((ir->eI == eiVV) || (ir->eI == eiVVAK)) && (ir->epc != epcMTTK) && (ir->etc == etcNOSEHOOVER));
1564 }
1565
1566 gmx_bool inputrecNphTrotter(const t_inputrec* ir)
1567 {
1568     return (((ir->eI == eiVV) || (ir->eI == eiVVAK)) && (ir->epc == epcMTTK) && (ir->etc != etcNOSEHOOVER));
1569 }
1570
1571 bool inputrecPbcXY2Walls(const t_inputrec* ir)
1572 {
1573     return (ir->pbcType == PbcType::XY && ir->nwall == 2);
1574 }
1575
1576 bool integratorHasConservedEnergyQuantity(const t_inputrec* ir)
1577 {
1578     if (!EI_MD(ir->eI))
1579     { // NOLINT bugprone-branch-clone
1580         // Energy minimization or stochastic integrator: no conservation
1581         return false;
1582     }
1583     else if (ir->etc == etcNO && ir->epc == epcNO)
1584     {
1585         // The total energy is conserved, no additional conserved quanitity
1586         return false;
1587     }
1588     else
1589     {
1590         // Shear stress with Parrinello-Rahman is not supported (tedious)
1591         bool shearWithPR =
1592                 ((ir->epc == epcPARRINELLORAHMAN || ir->epc == epcMTTK)
1593                  && (ir->ref_p[YY][XX] != 0 || ir->ref_p[ZZ][XX] != 0 || ir->ref_p[ZZ][YY] != 0));
1594
1595         return !ETC_ANDERSEN(ir->etc) && !shearWithPR;
1596     }
1597 }
1598
1599 bool integratorHasReferenceTemperature(const t_inputrec* ir)
1600 {
1601     return ((ir->etc != etcNO) || EI_SD(ir->eI) || (ir->eI == eiBD) || EI_TPI(ir->eI));
1602 }
1603
1604 int inputrec2nboundeddim(const t_inputrec* ir)
1605 {
1606     if (inputrecPbcXY2Walls(ir))
1607     {
1608         return 3;
1609     }
1610     else
1611     {
1612         return numPbcDimensions(ir->pbcType);
1613     }
1614 }
1615
1616 int ndof_com(const t_inputrec* ir)
1617 {
1618     int n = 0;
1619
1620     switch (ir->pbcType)
1621     {
1622         case PbcType::Xyz:
1623         case PbcType::No: n = 3; break;
1624         case PbcType::XY: n = (ir->nwall == 0 ? 3 : 2); break;
1625         case PbcType::Screw: n = 1; break;
1626         default: gmx_incons("Unknown pbc in calc_nrdf");
1627     }
1628
1629     return n;
1630 }
1631
1632 real maxReferenceTemperature(const t_inputrec& ir)
1633 {
1634     if (EI_ENERGY_MINIMIZATION(ir.eI) || ir.eI == eiNM)
1635     {
1636         return 0;
1637     }
1638
1639     if (EI_MD(ir.eI) && ir.etc == etcNO)
1640     {
1641         return -1;
1642     }
1643
1644     /* SD and BD also use ref_t and tau_t for setting the reference temperature.
1645      * TPI can be treated as MD, since it needs an ensemble temperature.
1646      */
1647
1648     real maxTemperature = 0;
1649     for (int i = 0; i < ir.opts.ngtc; i++)
1650     {
1651         if (ir.opts.tau_t[i] >= 0)
1652         {
1653             maxTemperature = std::max(maxTemperature, ir.opts.ref_t[i]);
1654         }
1655     }
1656
1657     return maxTemperature;
1658 }
1659
1660 bool haveEwaldSurfaceContribution(const t_inputrec& ir)
1661 {
1662     return EEL_PME_EWALD(ir.coulombtype) && (ir.ewald_geometry == eewg3DC || ir.epsilon_surface != 0);
1663 }
1664
1665 bool haveFreeEnergyType(const t_inputrec& ir, const int fepType)
1666 {
1667     for (int i = 0; i < ir.fepvals->n_lambda; i++)
1668     {
1669         if (ir.fepvals->all_lambda[fepType][i] > 0)
1670         {
1671             return true;
1672         }
1673     }
1674     return false;
1675 }