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