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