2 * This file is part of the GROMACS molecular simulation package.
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,2019, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
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.
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.
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.
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.
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.
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"
64 //! Macro to select a bool name
65 #define EBOOL(e) gmx::boolToString(e)
67 /* The minimum number of integration steps required for reasonably accurate
68 * integration of first and second order coupling algorithms.
70 const int nstmin_berendsen_tcouple = 5;
71 const int nstmin_berendsen_pcouple = 10;
72 const int nstmin_harmonic = 20;
74 t_inputrec::t_inputrec()
76 // TODO When this memset is removed, remove the suppression of
77 // gcc -Wno-class-memaccess in a CMakeLists.txt file.
78 std::memset(this, 0, sizeof(*this)); // NOLINT(bugprone-undefined-memory-manipulation)
80 snew(expandedvals, 1);
84 t_inputrec::~t_inputrec()
89 static int nst_wanted(const t_inputrec* ir)
101 int ir_optimal_nstcalcenergy(const t_inputrec* ir)
103 return nst_wanted(ir);
106 int tcouple_min_integration_steps(int etc)
112 case etcNO: n = 0; break;
114 case etcYES: n = nstmin_berendsen_tcouple; break;
116 /* V-rescale supports instantaneous rescaling */
119 case etcNOSEHOOVER: n = nstmin_harmonic; break;
121 case etcANDERSENMASSIVE: n = 1; break;
122 default: gmx_incons("Unknown etc value");
128 int ir_optimal_nsttcouple(const t_inputrec* ir)
130 int nmin, nwanted, n;
134 nmin = tcouple_min_integration_steps(ir->etc);
136 nwanted = nst_wanted(ir);
139 if (ir->etc != etcNO)
141 for (g = 0; g < ir->opts.ngtc; g++)
143 if (ir->opts.tau_t[g] > 0)
145 tau_min = std::min(tau_min, ir->opts.tau_t[g]);
150 if (nmin == 0 || ir->delta_t * nwanted <= tau_min)
156 n = static_cast<int>(tau_min / (ir->delta_t * nmin) + 0.001);
161 while (nwanted % n != 0)
170 int pcouple_min_integration_steps(int epc)
176 case epcNO: n = 0; break;
178 case epcISOTROPIC: n = nstmin_berendsen_pcouple; break;
179 case epcPARRINELLORAHMAN:
180 case epcMTTK: n = nstmin_harmonic; break;
181 default: gmx_incons("Unknown epc value");
187 int ir_optimal_nstpcouple(const t_inputrec* ir)
189 int nmin, nwanted, n;
191 nmin = pcouple_min_integration_steps(ir->epc);
193 nwanted = nst_wanted(ir);
195 if (nmin == 0 || ir->delta_t * nwanted <= ir->tau_p)
201 n = static_cast<int>(ir->tau_p / (ir->delta_t * nmin) + 0.001);
206 while (nwanted % n != 0)
215 gmx_bool ir_coulomb_switched(const t_inputrec* ir)
217 return (ir->coulombtype == eelSWITCH || ir->coulombtype == eelSHIFT || ir->coulombtype == eelENCADSHIFT
218 || ir->coulombtype == eelPMESWITCH || ir->coulombtype == eelPMEUSERSWITCH
219 || ir->coulomb_modifier == eintmodPOTSWITCH || ir->coulomb_modifier == eintmodFORCESWITCH);
222 gmx_bool ir_coulomb_is_zero_at_cutoff(const t_inputrec* ir)
224 return (ir->cutoff_scheme == ecutsVERLET || ir_coulomb_switched(ir)
225 || ir->coulomb_modifier != eintmodNONE || ir->coulombtype == eelRF_ZERO);
228 gmx_bool ir_coulomb_might_be_zero_at_cutoff(const t_inputrec* ir)
230 return (ir_coulomb_is_zero_at_cutoff(ir) || ir->coulombtype == eelUSER || ir->coulombtype == eelPMEUSER);
233 gmx_bool ir_vdw_switched(const t_inputrec* ir)
235 return (ir->vdwtype == evdwSWITCH || ir->vdwtype == evdwSHIFT || ir->vdwtype == evdwENCADSHIFT
236 || ir->vdw_modifier == eintmodPOTSWITCH || ir->vdw_modifier == eintmodFORCESWITCH);
239 gmx_bool ir_vdw_is_zero_at_cutoff(const t_inputrec* ir)
241 return (ir->cutoff_scheme == ecutsVERLET || ir_vdw_switched(ir) || ir->vdw_modifier != eintmodNONE);
244 gmx_bool ir_vdw_might_be_zero_at_cutoff(const t_inputrec* ir)
246 return (ir_vdw_is_zero_at_cutoff(ir) || ir->vdwtype == evdwUSER);
249 static void done_pull_group(t_pull_group* pgrp)
258 static void done_pull_params(pull_params_t* pull)
262 for (i = 0; i < pull->ngroup + 1; i++)
264 done_pull_group(pull->group);
271 static void done_lambdas(t_lambda* fep)
273 if (fep->n_lambda > 0)
275 for (int i = 0; i < efptNR; i++)
277 sfree(fep->all_lambda[i]);
280 sfree(fep->all_lambda);
283 void done_inputrec(t_inputrec* ir)
285 sfree(ir->opts.nrdf);
286 sfree(ir->opts.ref_t);
287 for (int i = 0; i < ir->opts.ngtc; i++)
289 sfree(ir->opts.anneal_time[i]);
290 sfree(ir->opts.anneal_temp[i]);
292 sfree(ir->opts.annealing);
293 sfree(ir->opts.anneal_npoints);
294 sfree(ir->opts.anneal_time);
295 sfree(ir->opts.anneal_temp);
296 sfree(ir->opts.tau_t);
298 sfree(ir->opts.nFreeze);
299 sfree(ir->opts.QMmethod);
300 sfree(ir->opts.QMbasis);
301 sfree(ir->opts.QMcharge);
302 sfree(ir->opts.QMmult);
304 sfree(ir->opts.CASorbitals);
305 sfree(ir->opts.CASelectrons);
306 sfree(ir->opts.SAon);
307 sfree(ir->opts.SAoff);
308 sfree(ir->opts.SAsteps);
309 sfree(ir->opts.egp_flags);
310 done_lambdas(ir->fepvals);
312 sfree(ir->expandedvals);
313 sfree(ir->simtempvals);
317 done_pull_params(ir->pull);
323 static void pr_qm_opts(FILE* fp, int indent, const char* title, const t_grpopts* opts)
325 fprintf(fp, "%s:\n", title);
327 pr_int(fp, indent, "ngQM", opts->ngQM);
330 pr_ivec(fp, indent, "QMmethod", opts->QMmethod, opts->ngQM, FALSE);
331 pr_ivec(fp, indent, "QMbasis", opts->QMbasis, opts->ngQM, FALSE);
332 pr_ivec(fp, indent, "QMcharge", opts->QMcharge, opts->ngQM, FALSE);
333 pr_ivec(fp, indent, "QMmult", opts->QMmult, opts->ngQM, FALSE);
334 pr_bvec(fp, indent, "SH", opts->bSH, opts->ngQM, FALSE);
335 pr_ivec(fp, indent, "CASorbitals", opts->CASorbitals, opts->ngQM, FALSE);
336 pr_ivec(fp, indent, "CASelectrons", opts->CASelectrons, opts->ngQM, FALSE);
337 pr_rvec(fp, indent, "SAon", opts->SAon, opts->ngQM, FALSE);
338 pr_rvec(fp, indent, "SAoff", opts->SAoff, opts->ngQM, FALSE);
339 pr_ivec(fp, indent, "SAsteps", opts->SAsteps, opts->ngQM, FALSE);
343 static void pr_grp_opts(FILE* out, int indent, const char* title, const t_grpopts* opts, gmx_bool bMDPformat)
349 fprintf(out, "%s:\n", title);
352 pr_indent(out, indent);
353 fprintf(out, "nrdf%s", bMDPformat ? " = " : ":");
354 for (i = 0; (i < opts->ngtc); i++)
356 fprintf(out, " %10g", opts->nrdf[i]);
360 pr_indent(out, indent);
361 fprintf(out, "ref-t%s", bMDPformat ? " = " : ":");
362 for (i = 0; (i < opts->ngtc); i++)
364 fprintf(out, " %10g", opts->ref_t[i]);
368 pr_indent(out, indent);
369 fprintf(out, "tau-t%s", bMDPformat ? " = " : ":");
370 for (i = 0; (i < opts->ngtc); i++)
372 fprintf(out, " %10g", opts->tau_t[i]);
376 /* Pretty-print the simulated annealing info */
377 fprintf(out, "annealing%s", bMDPformat ? " = " : ":");
378 for (i = 0; (i < opts->ngtc); i++)
380 fprintf(out, " %10s", EANNEAL(opts->annealing[i]));
384 fprintf(out, "annealing-npoints%s", bMDPformat ? " = " : ":");
385 for (i = 0; (i < opts->ngtc); i++)
387 fprintf(out, " %10d", opts->anneal_npoints[i]);
391 for (i = 0; (i < opts->ngtc); i++)
393 if (opts->anneal_npoints[i] > 0)
395 fprintf(out, "annealing-time [%d]:\t", i);
396 for (j = 0; (j < opts->anneal_npoints[i]); j++)
398 fprintf(out, " %10.1f", opts->anneal_time[i][j]);
401 fprintf(out, "annealing-temp [%d]:\t", i);
402 for (j = 0; (j < opts->anneal_npoints[i]); j++)
404 fprintf(out, " %10.1f", opts->anneal_temp[i][j]);
410 pr_indent(out, indent);
411 fprintf(out, "acc:\t");
412 for (i = 0; (i < opts->ngacc); i++)
414 for (m = 0; (m < DIM); m++)
416 fprintf(out, " %10g", opts->acc[i][m]);
421 pr_indent(out, indent);
422 fprintf(out, "nfreeze:");
423 for (i = 0; (i < opts->ngfrz); i++)
425 for (m = 0; (m < DIM); m++)
427 fprintf(out, " %10s", opts->nFreeze[i][m] ? "Y" : "N");
433 for (i = 0; (i < opts->ngener); i++)
435 pr_indent(out, indent);
436 fprintf(out, "energygrp-flags[%3d]:", i);
437 for (m = 0; (m < opts->ngener); m++)
439 fprintf(out, " %d", opts->egp_flags[opts->ngener * i + m]);
447 static void pr_matrix(FILE* fp, int indent, const char* title, const rvec* m, gmx_bool bMDPformat)
451 fprintf(fp, "%-10s = %g %g %g %g %g %g\n", title, m[XX][XX], m[YY][YY], m[ZZ][ZZ],
452 m[XX][YY], m[XX][ZZ], m[YY][ZZ]);
456 pr_rvecs(fp, indent, title, m, DIM);
460 #define PS(t, s) pr_str(fp, indent, t, s)
461 #define PI(t, s) pr_int(fp, indent, t, s)
462 #define PSTEP(t, s) pr_int64(fp, indent, t, s)
463 #define PR(t, s) pr_real(fp, indent, t, s)
464 #define PD(t, s) pr_double(fp, indent, t, s)
466 static void pr_pull_group(FILE* fp, int indent, int g, const t_pull_group* pgrp)
468 pr_indent(fp, indent);
469 fprintf(fp, "pull-group %d:\n", g);
471 pr_ivec_block(fp, indent, "atom", pgrp->ind, pgrp->nat, TRUE);
472 pr_rvec(fp, indent, "weight", pgrp->weight, pgrp->nweight, TRUE);
473 PI("pbcatom", pgrp->pbcatom);
476 static void pr_pull_coord(FILE* fp, int indent, int c, const t_pull_coord* pcrd)
480 pr_indent(fp, indent);
481 fprintf(fp, "pull-coord %d:\n", c);
482 PS("type", EPULLTYPE(pcrd->eType));
483 if (pcrd->eType == epullEXTERNAL)
485 PS("potential-provider", pcrd->externalPotentialProvider);
487 PS("geometry", EPULLGEOM(pcrd->eGeom));
488 for (g = 0; g < pcrd->ngroup; g++)
492 sprintf(buf, "group[%d]", g);
493 PI(buf, pcrd->group[g]);
495 pr_ivec(fp, indent, "dim", pcrd->dim, DIM, TRUE);
496 pr_rvec(fp, indent, "origin", pcrd->origin, DIM, TRUE);
497 pr_rvec(fp, indent, "vec", pcrd->vec, DIM, TRUE);
498 PS("start", EBOOL(pcrd->bStart));
499 PR("init", pcrd->init);
500 PR("rate", pcrd->rate);
505 static void pr_simtempvals(FILE* fp, int indent, const t_simtemp* simtemp, int n_lambda)
507 PS("simulated-tempering-scaling", ESIMTEMP(simtemp->eSimTempScale));
508 PR("sim-temp-low", simtemp->simtemp_low);
509 PR("sim-temp-high", simtemp->simtemp_high);
510 pr_rvec(fp, indent, "simulated tempering temperatures", simtemp->temperatures, n_lambda, TRUE);
513 static void pr_expandedvals(FILE* fp, int indent, const t_expanded* expand, int n_lambda)
516 PI("nstexpanded", expand->nstexpanded);
517 PS("lmc-stats", elamstats_names[expand->elamstats]);
518 PS("lmc-move", elmcmove_names[expand->elmcmove]);
519 PS("lmc-weights-equil", elmceq_names[expand->elmceq]);
520 if (expand->elmceq == elmceqNUMATLAM)
522 PI("weight-equil-number-all-lambda", expand->equil_n_at_lam);
524 if (expand->elmceq == elmceqSAMPLES)
526 PI("weight-equil-number-samples", expand->equil_samples);
528 if (expand->elmceq == elmceqSTEPS)
530 PI("weight-equil-number-steps", expand->equil_steps);
532 if (expand->elmceq == elmceqWLDELTA)
534 PR("weight-equil-wl-delta", expand->equil_wl_delta);
536 if (expand->elmceq == elmceqRATIO)
538 PR("weight-equil-count-ratio", expand->equil_ratio);
540 PI("lmc-seed", expand->lmc_seed);
541 PR("mc-temperature", expand->mc_temp);
542 PI("lmc-repeats", expand->lmc_repeats);
543 PI("lmc-gibbsdelta", expand->gibbsdeltalam);
544 PI("lmc-forced-nstart", expand->lmc_forced_nstart);
545 PS("symmetrized-transition-matrix", EBOOL(expand->bSymmetrizedTMatrix));
546 PI("nst-transition-matrix", expand->nstTij);
547 PI("mininum-var-min", expand->minvarmin); /*default is reasonable */
548 PI("weight-c-range", expand->c_range); /* default is just C=0 */
549 PR("wl-scale", expand->wl_scale);
550 PR("wl-ratio", expand->wl_ratio);
551 PR("init-wl-delta", expand->init_wl_delta);
552 PS("wl-oneovert", EBOOL(expand->bWLoneovert));
554 pr_indent(fp, indent);
555 pr_rvec(fp, indent, "init-lambda-weights", expand->init_lambda_weights, n_lambda, TRUE);
556 PS("init-weights", EBOOL(expand->bInit_weights));
559 static void pr_fepvals(FILE* fp, int indent, const t_lambda* fep, gmx_bool bMDPformat)
563 PR("init-lambda", fep->init_lambda);
564 PI("init-lambda-state", fep->init_fep_state);
565 PR("delta-lambda", fep->delta_lambda);
566 PI("nstdhdl", fep->nstdhdl);
570 PI("n-lambdas", fep->n_lambda);
572 if (fep->n_lambda > 0)
574 pr_indent(fp, indent);
575 fprintf(fp, "separate-dvdl%s\n", bMDPformat ? " = " : ":");
576 for (i = 0; i < efptNR; i++)
578 fprintf(fp, "%18s = ", efpt_names[i]);
579 if (fep->separate_dvdl[i])
581 fprintf(fp, " TRUE");
585 fprintf(fp, " FALSE");
589 fprintf(fp, "all-lambdas%s\n", bMDPformat ? " = " : ":");
590 for (i = 0; i < efptNR; i++)
592 fprintf(fp, "%18s = ", efpt_names[i]);
593 for (j = 0; j < fep->n_lambda; j++)
595 fprintf(fp, " %10g", fep->all_lambda[i][j]);
600 PI("calc-lambda-neighbors", fep->lambda_neighbors);
601 PS("dhdl-print-energy", edHdLPrintEnergy_names[fep->edHdLPrintEnergy]);
602 PR("sc-alpha", fep->sc_alpha);
603 PI("sc-power", fep->sc_power);
604 PR("sc-r-power", fep->sc_r_power);
605 PR("sc-sigma", fep->sc_sigma);
606 PR("sc-sigma-min", fep->sc_sigma_min);
607 PS("sc-coul", EBOOL(fep->bScCoul));
608 PI("dh-hist-size", fep->dh_hist_size);
609 PD("dh-hist-spacing", fep->dh_hist_spacing);
610 PS("separate-dhdl-file", SEPDHDLFILETYPE(fep->separate_dhdl_file));
611 PS("dhdl-derivatives", DHDLDERIVATIVESTYPE(fep->dhdl_derivatives));
614 static void pr_pull(FILE* fp, int indent, const pull_params_t* pull)
618 PR("pull-cylinder-r", pull->cylinder_r);
619 PR("pull-constr-tol", pull->constr_tol);
620 PS("pull-print-COM", EBOOL(pull->bPrintCOM));
621 PS("pull-print-ref-value", EBOOL(pull->bPrintRefValue));
622 PS("pull-print-components", EBOOL(pull->bPrintComp));
623 PI("pull-nstxout", pull->nstxout);
624 PI("pull-nstfout", pull->nstfout);
625 PS("pull-pbc-ref-prev-step-com", EBOOL(pull->bSetPbcRefToPrevStepCOM));
626 PS("pull-xout-average", EBOOL(pull->bXOutAverage));
627 PS("pull-fout-average", EBOOL(pull->bFOutAverage));
628 PI("pull-ngroups", pull->ngroup);
629 for (g = 0; g < pull->ngroup; g++)
631 pr_pull_group(fp, indent, g, &pull->group[g]);
633 PI("pull-ncoords", pull->ncoord);
634 for (g = 0; g < pull->ncoord; g++)
636 pr_pull_coord(fp, indent, g, &pull->coord[g]);
640 static void pr_awh_bias_dim(FILE* fp, int indent, gmx::AwhDimParams* awhDimParams, const char* prefix)
642 pr_indent(fp, indent);
644 fprintf(fp, "%s:\n", prefix);
645 PS("coord-provider", EAWHCOORDPROVIDER(awhDimParams->eCoordProvider));
646 PI("coord-index", awhDimParams->coordIndex + 1);
647 PR("start", awhDimParams->origin);
648 PR("end", awhDimParams->end);
649 PR("period", awhDimParams->period);
650 PR("force-constant", awhDimParams->forceConstant);
651 PR("diffusion", awhDimParams->diffusion);
652 PR("start", awhDimParams->origin);
653 PR("end", awhDimParams->end);
654 PR("cover-diameter", awhDimParams->coverDiameter);
657 static void pr_awh_bias(FILE* fp, int indent, gmx::AwhBiasParams* awhBiasParams, const char* prefix)
661 sprintf(opt, "%s-error-init", prefix);
662 PR(opt, awhBiasParams->errorInitial);
663 sprintf(opt, "%s-growth", prefix);
664 PS(opt, EAWHGROWTH(awhBiasParams->eGrowth));
665 sprintf(opt, "%s-target", prefix);
666 PS(opt, EAWHTARGET(awhBiasParams->eTarget));
667 sprintf(opt, "%s-target-beta-scalng", prefix);
668 PR(opt, awhBiasParams->targetBetaScaling);
669 sprintf(opt, "%s-target-cutoff", prefix);
670 PR(opt, awhBiasParams->targetCutoff);
671 sprintf(opt, "%s-user-data", prefix);
672 PS(opt, EBOOL(awhBiasParams->bUserData));
673 sprintf(opt, "%s-share-group", prefix);
674 PI(opt, awhBiasParams->shareGroup);
675 sprintf(opt, "%s-equilibrate-histogram", prefix);
676 PS(opt, EBOOL(awhBiasParams->equilibrateHistogram));
677 sprintf(opt, "%s-ndim", prefix);
678 PI(opt, awhBiasParams->ndim);
680 for (int d = 0; d < awhBiasParams->ndim; d++)
682 char prefixdim[STRLEN];
683 sprintf(prefixdim, "%s-dim%d", prefix, d + 1);
684 pr_awh_bias_dim(fp, indent, &awhBiasParams->dimParams[d], prefixdim);
688 static void pr_awh(FILE* fp, int indent, gmx::AwhParams* awhParams)
690 PS("awh-potential", EAWHPOTENTIAL(awhParams->ePotential));
691 PI("awh-seed", awhParams->seed);
692 PI("awh-nstout", awhParams->nstOut);
693 PI("awh-nstsample", awhParams->nstSampleCoord);
694 PI("awh-nsamples-update", awhParams->numSamplesUpdateFreeEnergy);
695 PS("awh-share-bias-multisim", EBOOL(awhParams->shareBiasMultisim));
696 PI("awh-nbias", awhParams->numBias);
698 for (int k = 0; k < awhParams->numBias; k++)
700 auto prefix = gmx::formatString("awh%d", k + 1);
701 pr_awh_bias(fp, indent, &awhParams->awhBiasParams[k], prefix.c_str());
705 static void pr_rotgrp(FILE* fp, int indent, int g, const t_rotgrp* rotg)
707 pr_indent(fp, indent);
708 fprintf(fp, "rot-group %d:\n", g);
710 PS("rot-type", EROTGEOM(rotg->eType));
711 PS("rot-massw", EBOOL(rotg->bMassW));
712 pr_ivec_block(fp, indent, "atom", rotg->ind, rotg->nat, TRUE);
713 pr_rvecs(fp, indent, "x-ref", rotg->x_ref, rotg->nat);
714 pr_rvec(fp, indent, "rot-vec", rotg->inputVec, DIM, TRUE);
715 pr_rvec(fp, indent, "rot-pivot", rotg->pivot, DIM, TRUE);
716 PR("rot-rate", rotg->rate);
717 PR("rot-k", rotg->k);
718 PR("rot-slab-dist", rotg->slab_dist);
719 PR("rot-min-gauss", rotg->min_gaussian);
720 PR("rot-eps", rotg->eps);
721 PS("rot-fit-method", EROTFIT(rotg->eFittype));
722 PI("rot-potfit-nstep", rotg->PotAngle_nstep);
723 PR("rot-potfit-step", rotg->PotAngle_step);
726 static void pr_rot(FILE* fp, int indent, const t_rot* rot)
730 PI("rot-nstrout", rot->nstrout);
731 PI("rot-nstsout", rot->nstsout);
732 PI("rot-ngroups", rot->ngrp);
733 for (g = 0; g < rot->ngrp; g++)
735 pr_rotgrp(fp, indent, g, &rot->grp[g]);
740 static void pr_swap(FILE* fp, int indent, const t_swapcoords* swap)
744 /* Enums for better readability of the code */
752 PI("swap-frequency", swap->nstswap);
754 /* The split groups that define the compartments */
755 for (int j = 0; j < 2; j++)
757 snprintf(str, STRLEN, "massw_split%d", j);
758 PS(str, EBOOL(swap->massw_split[j]));
759 snprintf(str, STRLEN, "split atoms group %d", j);
760 pr_ivec_block(fp, indent, str, swap->grp[j].ind, swap->grp[j].nat, TRUE);
763 /* The solvent group */
764 snprintf(str, STRLEN, "solvent group %s", swap->grp[eGrpSolvent].molname);
765 pr_ivec_block(fp, indent, str, swap->grp[eGrpSolvent].ind, swap->grp[eGrpSolvent].nat, TRUE);
767 /* Now print the indices for all the ion groups: */
768 for (int ig = eSwapFixedGrpNR; ig < swap->ngrp; ig++)
770 snprintf(str, STRLEN, "ion group %s", swap->grp[ig].molname);
771 pr_ivec_block(fp, indent, str, swap->grp[ig].ind, swap->grp[ig].nat, TRUE);
774 PR("cyl0-r", swap->cyl0r);
775 PR("cyl0-up", swap->cyl0u);
776 PR("cyl0-down", swap->cyl0l);
777 PR("cyl1-r", swap->cyl1r);
778 PR("cyl1-up", swap->cyl1u);
779 PR("cyl1-down", swap->cyl1l);
780 PI("coupl-steps", swap->nAverage);
782 /* Print the requested ion counts for both compartments */
783 for (int ic = eCompA; ic <= eCompB; ic++)
785 for (int ig = eSwapFixedGrpNR; ig < swap->ngrp; ig++)
787 snprintf(str, STRLEN, "%s-in-%c", swap->grp[ig].molname, 'A' + ic);
788 PI(str, swap->grp[ig].nmolReq[ic]);
792 PR("threshold", swap->threshold);
793 PR("bulk-offsetA", swap->bulkOffset[eCompA]);
794 PR("bulk-offsetB", swap->bulkOffset[eCompB]);
798 static void pr_imd(FILE* fp, int indent, const t_IMD* imd)
800 PI("IMD-atoms", imd->nat);
801 pr_ivec_block(fp, indent, "atom", imd->ind, imd->nat, TRUE);
805 void pr_inputrec(FILE* fp, int indent, const char* title, const t_inputrec* ir, gmx_bool bMDPformat)
807 const char* infbuf = "inf";
809 if (available(fp, ir, indent, title))
813 indent = pr_title(fp, indent, title);
815 /* Try to make this list appear in the same order as the
816 * options are written in the default mdout.mdp, and with
817 * the same user-exposed names to facilitate debugging.
819 PS("integrator", EI(ir->eI));
820 PR("tinit", ir->init_t);
821 PR("dt", ir->delta_t);
822 PSTEP("nsteps", ir->nsteps);
823 PSTEP("init-step", ir->init_step);
824 PI("simulation-part", ir->simulation_part);
825 PS("comm-mode", ECOM(ir->comm_mode));
826 PI("nstcomm", ir->nstcomm);
828 /* Langevin dynamics */
829 PR("bd-fric", ir->bd_fric);
830 PSTEP("ld-seed", ir->ld_seed);
832 /* Energy minimization */
833 PR("emtol", ir->em_tol);
834 PR("emstep", ir->em_stepsize);
835 PI("niter", ir->niter);
836 PR("fcstep", ir->fc_stepsize);
837 PI("nstcgsteep", ir->nstcgsteep);
838 PI("nbfgscorr", ir->nbfgscorr);
840 /* Test particle insertion */
841 PR("rtpi", ir->rtpi);
844 PI("nstxout", ir->nstxout);
845 PI("nstvout", ir->nstvout);
846 PI("nstfout", ir->nstfout);
847 PI("nstlog", ir->nstlog);
848 PI("nstcalcenergy", ir->nstcalcenergy);
849 PI("nstenergy", ir->nstenergy);
850 PI("nstxout-compressed", ir->nstxout_compressed);
851 PR("compressed-x-precision", ir->x_compression_precision);
853 /* Neighborsearching parameters */
854 PS("cutoff-scheme", ECUTSCHEME(ir->cutoff_scheme));
855 PI("nstlist", ir->nstlist);
856 PS("pbc", epbc_names[ir->ePBC]);
857 PS("periodic-molecules", EBOOL(ir->bPeriodicMols));
858 PR("verlet-buffer-tolerance", ir->verletbuf_tol);
859 PR("rlist", ir->rlist);
861 /* Options for electrostatics and VdW */
862 PS("coulombtype", EELTYPE(ir->coulombtype));
863 PS("coulomb-modifier", INTMODIFIER(ir->coulomb_modifier));
864 PR("rcoulomb-switch", ir->rcoulomb_switch);
865 PR("rcoulomb", ir->rcoulomb);
866 if (ir->epsilon_r != 0)
868 PR("epsilon-r", ir->epsilon_r);
872 PS("epsilon-r", infbuf);
874 if (ir->epsilon_rf != 0)
876 PR("epsilon-rf", ir->epsilon_rf);
880 PS("epsilon-rf", infbuf);
882 PS("vdw-type", EVDWTYPE(ir->vdwtype));
883 PS("vdw-modifier", INTMODIFIER(ir->vdw_modifier));
884 PR("rvdw-switch", ir->rvdw_switch);
885 PR("rvdw", ir->rvdw);
886 PS("DispCorr", EDISPCORR(ir->eDispCorr));
887 PR("table-extension", ir->tabext);
889 PR("fourierspacing", ir->fourier_spacing);
890 PI("fourier-nx", ir->nkx);
891 PI("fourier-ny", ir->nky);
892 PI("fourier-nz", ir->nkz);
893 PI("pme-order", ir->pme_order);
894 PR("ewald-rtol", ir->ewald_rtol);
895 PR("ewald-rtol-lj", ir->ewald_rtol_lj);
896 PS("lj-pme-comb-rule", ELJPMECOMBNAMES(ir->ljpme_combination_rule));
897 PR("ewald-geometry", ir->ewald_geometry);
898 PR("epsilon-surface", ir->epsilon_surface);
900 /* Options for weak coupling algorithms */
901 PS("tcoupl", ETCOUPLTYPE(ir->etc));
902 PI("nsttcouple", ir->nsttcouple);
903 PI("nh-chain-length", ir->opts.nhchainlength);
904 PS("print-nose-hoover-chain-variables", EBOOL(ir->bPrintNHChains));
906 PS("pcoupl", EPCOUPLTYPE(ir->epc));
907 PS("pcoupltype", EPCOUPLTYPETYPE(ir->epct));
908 PI("nstpcouple", ir->nstpcouple);
909 PR("tau-p", ir->tau_p);
910 pr_matrix(fp, indent, "compressibility", ir->compress, bMDPformat);
911 pr_matrix(fp, indent, "ref-p", ir->ref_p, bMDPformat);
912 PS("refcoord-scaling", EREFSCALINGTYPE(ir->refcoord_scaling));
916 fprintf(fp, "posres-com = %g %g %g\n", ir->posres_com[XX], ir->posres_com[YY],
918 fprintf(fp, "posres-comB = %g %g %g\n", ir->posres_comB[XX], ir->posres_comB[YY],
919 ir->posres_comB[ZZ]);
923 pr_rvec(fp, indent, "posres-com", ir->posres_com, DIM, TRUE);
924 pr_rvec(fp, indent, "posres-comB", ir->posres_comB, DIM, TRUE);
928 PS("QMMM", EBOOL(ir->bQMMM));
929 PI("QMconstraints", ir->QMconstraints);
930 PI("QMMMscheme", ir->QMMMscheme);
931 PR("MMChargeScaleFactor", ir->scalefactor);
932 pr_qm_opts(fp, indent, "qm-opts", &(ir->opts));
934 /* CONSTRAINT OPTIONS */
935 PS("constraint-algorithm", ECONSTRTYPE(ir->eConstrAlg));
936 PS("continuation", EBOOL(ir->bContinuation));
938 PS("Shake-SOR", EBOOL(ir->bShakeSOR));
939 PR("shake-tol", ir->shake_tol);
940 PI("lincs-order", ir->nProjOrder);
941 PI("lincs-iter", ir->nLincsIter);
942 PR("lincs-warnangle", ir->LincsWarnAngle);
945 PI("nwall", ir->nwall);
946 PS("wall-type", EWALLTYPE(ir->wall_type));
947 PR("wall-r-linpot", ir->wall_r_linpot);
949 PI("wall-atomtype[0]", ir->wall_atomtype[0]);
950 PI("wall-atomtype[1]", ir->wall_atomtype[1]);
952 PR("wall-density[0]", ir->wall_density[0]);
953 PR("wall-density[1]", ir->wall_density[1]);
954 PR("wall-ewald-zfac", ir->wall_ewald_zfac);
957 PS("pull", EBOOL(ir->bPull));
960 pr_pull(fp, indent, ir->pull);
964 PS("awh", EBOOL(ir->bDoAwh));
967 pr_awh(fp, indent, ir->awhParams);
970 /* ENFORCED ROTATION */
971 PS("rotation", EBOOL(ir->bRot));
974 pr_rot(fp, indent, ir->rot);
978 PS("interactiveMD", EBOOL(ir->bIMD));
981 pr_imd(fp, indent, ir->imd);
984 /* NMR refinement stuff */
985 PS("disre", EDISRETYPE(ir->eDisre));
986 PS("disre-weighting", EDISREWEIGHTING(ir->eDisreWeighting));
987 PS("disre-mixed", EBOOL(ir->bDisreMixed));
988 PR("dr-fc", ir->dr_fc);
989 PR("dr-tau", ir->dr_tau);
990 PR("nstdisreout", ir->nstdisreout);
992 PR("orire-fc", ir->orires_fc);
993 PR("orire-tau", ir->orires_tau);
994 PR("nstorireout", ir->nstorireout);
996 /* FREE ENERGY VARIABLES */
997 PS("free-energy", EFEPTYPE(ir->efep));
998 if (ir->efep != efepNO || ir->bSimTemp)
1000 pr_fepvals(fp, indent, ir->fepvals, bMDPformat);
1004 pr_expandedvals(fp, indent, ir->expandedvals, ir->fepvals->n_lambda);
1007 /* NON-equilibrium MD stuff */
1008 PR("cos-acceleration", ir->cos_accel);
1009 pr_matrix(fp, indent, "deform", ir->deform, bMDPformat);
1011 /* SIMULATED TEMPERING */
1012 PS("simulated-tempering", EBOOL(ir->bSimTemp));
1015 pr_simtempvals(fp, indent, ir->simtempvals, ir->fepvals->n_lambda);
1018 /* ION/WATER SWAPPING FOR COMPUTATIONAL ELECTROPHYSIOLOGY */
1019 PS("swapcoords", ESWAPTYPE(ir->eSwapCoords));
1020 if (ir->eSwapCoords != eswapNO)
1022 pr_swap(fp, indent, ir->swap);
1025 /* USER-DEFINED THINGIES */
1026 PI("userint1", ir->userint1);
1027 PI("userint2", ir->userint2);
1028 PI("userint3", ir->userint3);
1029 PI("userint4", ir->userint4);
1030 PR("userreal1", ir->userreal1);
1031 PR("userreal2", ir->userreal2);
1032 PR("userreal3", ir->userreal3);
1033 PR("userreal4", ir->userreal4);
1037 gmx::TextWriter writer(fp);
1038 writer.wrapperSettings().setIndent(indent);
1039 gmx::dumpKeyValueTree(&writer, *ir->params);
1042 pr_grp_opts(fp, indent, "grpopts", &(ir->opts), bMDPformat);
1049 static void cmp_grpopts(FILE* fp, const t_grpopts* opt1, const t_grpopts* opt2, real ftol, real abstol)
1052 char buf1[256], buf2[256];
1054 cmp_int(fp, "inputrec->grpopts.ngtc", -1, opt1->ngtc, opt2->ngtc);
1055 cmp_int(fp, "inputrec->grpopts.ngacc", -1, opt1->ngacc, opt2->ngacc);
1056 cmp_int(fp, "inputrec->grpopts.ngfrz", -1, opt1->ngfrz, opt2->ngfrz);
1057 cmp_int(fp, "inputrec->grpopts.ngener", -1, opt1->ngener, opt2->ngener);
1058 for (i = 0; (i < std::min(opt1->ngtc, opt2->ngtc)); i++)
1060 cmp_real(fp, "inputrec->grpopts.nrdf", i, opt1->nrdf[i], opt2->nrdf[i], ftol, abstol);
1061 cmp_real(fp, "inputrec->grpopts.ref_t", i, opt1->ref_t[i], opt2->ref_t[i], ftol, abstol);
1062 cmp_real(fp, "inputrec->grpopts.tau_t", i, opt1->tau_t[i], opt2->tau_t[i], ftol, abstol);
1063 cmp_int(fp, "inputrec->grpopts.annealing", i, opt1->annealing[i], opt2->annealing[i]);
1064 cmp_int(fp, "inputrec->grpopts.anneal_npoints", i, opt1->anneal_npoints[i],
1065 opt2->anneal_npoints[i]);
1066 if (opt1->anneal_npoints[i] == opt2->anneal_npoints[i])
1068 sprintf(buf1, "inputrec->grpopts.anneal_time[%d]", i);
1069 sprintf(buf2, "inputrec->grpopts.anneal_temp[%d]", i);
1070 for (j = 0; j < opt1->anneal_npoints[i]; j++)
1072 cmp_real(fp, buf1, j, opt1->anneal_time[i][j], opt2->anneal_time[i][j], ftol, abstol);
1073 cmp_real(fp, buf2, j, opt1->anneal_temp[i][j], opt2->anneal_temp[i][j], ftol, abstol);
1077 if (opt1->ngener == opt2->ngener)
1079 for (i = 0; i < opt1->ngener; i++)
1081 for (j = i; j < opt1->ngener; j++)
1083 sprintf(buf1, "inputrec->grpopts.egp_flags[%d]", i);
1084 cmp_int(fp, buf1, j, opt1->egp_flags[opt1->ngener * i + j],
1085 opt2->egp_flags[opt1->ngener * i + j]);
1089 for (i = 0; (i < std::min(opt1->ngacc, opt2->ngacc)); i++)
1091 cmp_rvec(fp, "inputrec->grpopts.acc", i, opt1->acc[i], opt2->acc[i], ftol, abstol);
1093 for (i = 0; (i < std::min(opt1->ngfrz, opt2->ngfrz)); i++)
1095 cmp_ivec(fp, "inputrec->grpopts.nFreeze", i, opt1->nFreeze[i], opt2->nFreeze[i]);
1099 static void cmp_pull(FILE* fp)
1102 "WARNING: Both files use COM pulling, but comparing of the pull struct is not "
1103 "implemented (yet). The pull parameters could be the same or different.\n");
1106 static void cmp_awhDimParams(FILE* fp,
1107 const gmx::AwhDimParams* dimp1,
1108 const gmx::AwhDimParams* dimp2,
1113 /* Note that we have double index here, but the compare functions only
1114 * support one index, so here we only print the dim index and not the bias.
1116 cmp_int(fp, "inputrec.awhParams->bias?->dim->coord_index", dimIndex, dimp1->coordIndex,
1118 cmp_double(fp, "inputrec->awhParams->bias?->dim->period", dimIndex, dimp1->period,
1119 dimp2->period, ftol, abstol);
1120 cmp_double(fp, "inputrec->awhParams->bias?->dim->diffusion", dimIndex, dimp1->diffusion,
1121 dimp2->diffusion, ftol, abstol);
1122 cmp_double(fp, "inputrec->awhParams->bias?->dim->origin", dimIndex, dimp1->origin,
1123 dimp2->origin, ftol, abstol);
1124 cmp_double(fp, "inputrec->awhParams->bias?->dim->end", dimIndex, dimp1->end, dimp2->end, ftol, abstol);
1125 cmp_double(fp, "inputrec->awhParams->bias?->dim->coord_value_init", dimIndex,
1126 dimp1->coordValueInit, dimp2->coordValueInit, ftol, abstol);
1127 cmp_double(fp, "inputrec->awhParams->bias?->dim->coverDiameter", dimIndex, dimp1->coverDiameter,
1128 dimp2->coverDiameter, ftol, abstol);
1131 static void cmp_awhBiasParams(FILE* fp,
1132 const gmx::AwhBiasParams* bias1,
1133 const gmx::AwhBiasParams* bias2,
1138 cmp_int(fp, "inputrec->awhParams->ndim", biasIndex, bias1->ndim, bias2->ndim);
1139 cmp_int(fp, "inputrec->awhParams->biaseTarget", biasIndex, bias1->eTarget, bias2->eTarget);
1140 cmp_double(fp, "inputrec->awhParams->biastargetBetaScaling", biasIndex,
1141 bias1->targetBetaScaling, bias2->targetBetaScaling, ftol, abstol);
1142 cmp_double(fp, "inputrec->awhParams->biastargetCutoff", biasIndex, bias1->targetCutoff,
1143 bias2->targetCutoff, ftol, abstol);
1144 cmp_int(fp, "inputrec->awhParams->biaseGrowth", biasIndex, bias1->eGrowth, bias2->eGrowth);
1145 cmp_bool(fp, "inputrec->awhParams->biasbUserData", biasIndex, bias1->bUserData != 0,
1146 bias2->bUserData != 0);
1147 cmp_double(fp, "inputrec->awhParams->biaserror_initial", biasIndex, bias1->errorInitial,
1148 bias2->errorInitial, ftol, abstol);
1149 cmp_int(fp, "inputrec->awhParams->biasShareGroup", biasIndex, bias1->shareGroup, bias2->shareGroup);
1151 for (int dim = 0; dim < std::min(bias1->ndim, bias2->ndim); dim++)
1153 cmp_awhDimParams(fp, &bias1->dimParams[dim], &bias2->dimParams[dim], dim, ftol, abstol);
1157 static void cmp_awhParams(FILE* fp, const gmx::AwhParams* awh1, const gmx::AwhParams* awh2, real ftol, real abstol)
1159 cmp_int(fp, "inputrec->awhParams->nbias", -1, awh1->numBias, awh2->numBias);
1160 cmp_int64(fp, "inputrec->awhParams->seed", awh1->seed, awh2->seed);
1161 cmp_int(fp, "inputrec->awhParams->nstout", -1, awh1->nstOut, awh2->nstOut);
1162 cmp_int(fp, "inputrec->awhParams->nstsample_coord", -1, awh1->nstSampleCoord, awh2->nstSampleCoord);
1163 cmp_int(fp, "inputrec->awhParams->nsamples_update_free_energy", -1,
1164 awh1->numSamplesUpdateFreeEnergy, awh2->numSamplesUpdateFreeEnergy);
1165 cmp_int(fp, "inputrec->awhParams->ePotential", -1, awh1->ePotential, awh2->ePotential);
1166 cmp_bool(fp, "inputrec->awhParams->shareBiasMultisim", -1, awh1->shareBiasMultisim,
1167 awh2->shareBiasMultisim);
1169 if (awh1->numBias == awh2->numBias)
1171 for (int bias = 0; bias < awh1->numBias; bias++)
1173 cmp_awhBiasParams(fp, &awh1->awhBiasParams[bias], &awh2->awhBiasParams[bias], bias, ftol, abstol);
1178 static void cmp_simtempvals(FILE* fp,
1179 const t_simtemp* simtemp1,
1180 const t_simtemp* simtemp2,
1186 cmp_int(fp, "inputrec->simtempvals->eSimTempScale", -1, simtemp1->eSimTempScale, simtemp2->eSimTempScale);
1187 cmp_real(fp, "inputrec->simtempvals->simtemp_high", -1, simtemp1->simtemp_high,
1188 simtemp2->simtemp_high, ftol, abstol);
1189 cmp_real(fp, "inputrec->simtempvals->simtemp_low", -1, simtemp1->simtemp_low,
1190 simtemp2->simtemp_low, ftol, abstol);
1191 for (i = 0; i < n_lambda; i++)
1193 cmp_real(fp, "inputrec->simtempvals->temperatures", -1, simtemp1->temperatures[i],
1194 simtemp2->temperatures[i], ftol, abstol);
1198 static void cmp_expandedvals(FILE* fp,
1199 const t_expanded* expand1,
1200 const t_expanded* expand2,
1207 cmp_bool(fp, "inputrec->fepvals->bInit_weights", -1, expand1->bInit_weights, expand2->bInit_weights);
1208 cmp_bool(fp, "inputrec->fepvals->bWLoneovert", -1, expand1->bWLoneovert, expand2->bWLoneovert);
1210 for (i = 0; i < n_lambda; i++)
1212 cmp_real(fp, "inputrec->expandedvals->init_lambda_weights", -1,
1213 expand1->init_lambda_weights[i], expand2->init_lambda_weights[i], ftol, abstol);
1216 cmp_int(fp, "inputrec->expandedvals->lambda-stats", -1, expand1->elamstats, expand2->elamstats);
1217 cmp_int(fp, "inputrec->expandedvals->lambda-mc-move", -1, expand1->elmcmove, expand2->elmcmove);
1218 cmp_int(fp, "inputrec->expandedvals->lmc-repeats", -1, expand1->lmc_repeats, expand2->lmc_repeats);
1219 cmp_int(fp, "inputrec->expandedvals->lmc-gibbsdelta", -1, expand1->gibbsdeltalam, expand2->gibbsdeltalam);
1220 cmp_int(fp, "inputrec->expandedvals->lmc-forced-nstart", -1, expand1->lmc_forced_nstart,
1221 expand2->lmc_forced_nstart);
1222 cmp_int(fp, "inputrec->expandedvals->lambda-weights-equil", -1, expand1->elmceq, expand2->elmceq);
1223 cmp_int(fp, "inputrec->expandedvals->,weight-equil-number-all-lambda", -1,
1224 expand1->equil_n_at_lam, expand2->equil_n_at_lam);
1225 cmp_int(fp, "inputrec->expandedvals->weight-equil-number-samples", -1, expand1->equil_samples,
1226 expand2->equil_samples);
1227 cmp_int(fp, "inputrec->expandedvals->weight-equil-number-steps", -1, expand1->equil_steps,
1228 expand2->equil_steps);
1229 cmp_real(fp, "inputrec->expandedvals->weight-equil-wl-delta", -1, expand1->equil_wl_delta,
1230 expand2->equil_wl_delta, ftol, abstol);
1231 cmp_real(fp, "inputrec->expandedvals->weight-equil-count-ratio", -1, expand1->equil_ratio,
1232 expand2->equil_ratio, ftol, abstol);
1233 cmp_bool(fp, "inputrec->expandedvals->symmetrized-transition-matrix", -1,
1234 expand1->bSymmetrizedTMatrix, expand2->bSymmetrizedTMatrix);
1235 cmp_int(fp, "inputrec->expandedvals->nstTij", -1, expand1->nstTij, expand2->nstTij);
1236 cmp_int(fp, "inputrec->expandedvals->mininum-var-min", -1, expand1->minvarmin,
1237 expand2->minvarmin); /*default is reasonable */
1238 cmp_int(fp, "inputrec->expandedvals->weight-c-range", -1, expand1->c_range, expand2->c_range); /* default is just C=0 */
1239 cmp_real(fp, "inputrec->expandedvals->wl-scale", -1, expand1->wl_scale, expand2->wl_scale, ftol, abstol);
1240 cmp_real(fp, "inputrec->expandedvals->init-wl-delta", -1, expand1->init_wl_delta,
1241 expand2->init_wl_delta, ftol, abstol);
1242 cmp_real(fp, "inputrec->expandedvals->wl-ratio", -1, expand1->wl_ratio, expand2->wl_ratio, ftol, abstol);
1243 cmp_int(fp, "inputrec->expandedvals->nstexpanded", -1, expand1->nstexpanded, expand2->nstexpanded);
1244 cmp_int(fp, "inputrec->expandedvals->lmc-seed", -1, expand1->lmc_seed, expand2->lmc_seed);
1245 cmp_real(fp, "inputrec->expandedvals->mc-temperature", -1, expand1->mc_temp, expand2->mc_temp,
1249 static void cmp_fepvals(FILE* fp, const t_lambda* fep1, const t_lambda* fep2, real ftol, real abstol)
1252 cmp_int(fp, "inputrec->nstdhdl", -1, fep1->nstdhdl, fep2->nstdhdl);
1253 cmp_double(fp, "inputrec->fepvals->init_fep_state", -1, fep1->init_fep_state,
1254 fep2->init_fep_state, ftol, abstol);
1255 cmp_double(fp, "inputrec->fepvals->delta_lambda", -1, fep1->delta_lambda, fep2->delta_lambda,
1257 cmp_int(fp, "inputrec->fepvals->n_lambda", -1, fep1->n_lambda, fep2->n_lambda);
1258 for (i = 0; i < efptNR; i++)
1260 for (j = 0; j < std::min(fep1->n_lambda, fep2->n_lambda); j++)
1262 cmp_double(fp, "inputrec->fepvals->all_lambda", -1, fep1->all_lambda[i][j],
1263 fep2->all_lambda[i][j], ftol, abstol);
1266 cmp_int(fp, "inputrec->fepvals->lambda_neighbors", 1, fep1->lambda_neighbors, fep2->lambda_neighbors);
1267 cmp_real(fp, "inputrec->fepvals->sc_alpha", -1, fep1->sc_alpha, fep2->sc_alpha, ftol, abstol);
1268 cmp_int(fp, "inputrec->fepvals->sc_power", -1, fep1->sc_power, fep2->sc_power);
1269 cmp_real(fp, "inputrec->fepvals->sc_r_power", -1, fep1->sc_r_power, fep2->sc_r_power, ftol, abstol);
1270 cmp_real(fp, "inputrec->fepvals->sc_sigma", -1, fep1->sc_sigma, fep2->sc_sigma, ftol, abstol);
1271 cmp_int(fp, "inputrec->fepvals->edHdLPrintEnergy", -1, fep1->edHdLPrintEnergy, fep1->edHdLPrintEnergy);
1272 cmp_bool(fp, "inputrec->fepvals->bScCoul", -1, fep1->bScCoul, fep1->bScCoul);
1273 cmp_int(fp, "inputrec->separate_dhdl_file", -1, fep1->separate_dhdl_file, fep2->separate_dhdl_file);
1274 cmp_int(fp, "inputrec->dhdl_derivatives", -1, fep1->dhdl_derivatives, fep2->dhdl_derivatives);
1275 cmp_int(fp, "inputrec->dh_hist_size", -1, fep1->dh_hist_size, fep2->dh_hist_size);
1276 cmp_double(fp, "inputrec->dh_hist_spacing", -1, fep1->dh_hist_spacing, fep2->dh_hist_spacing,
1280 void cmp_inputrec(FILE* fp, const t_inputrec* ir1, const t_inputrec* ir2, real ftol, real abstol)
1282 fprintf(fp, "comparing inputrec\n");
1284 /* gcc 2.96 doesnt like these defines at all, but issues a huge list
1285 * of warnings. Maybe it will change in future versions, but for the
1286 * moment I've spelled them out instead. /EL 000820
1287 * #define CIB(s) cmp_int(fp,"inputrec->"#s,0,ir1->##s,ir2->##s)
1288 * #define CII(s) cmp_int(fp,"inputrec->"#s,0,ir1->##s,ir2->##s)
1289 * #define CIR(s) cmp_real(fp,"inputrec->"#s,0,ir1->##s,ir2->##s,ftol)
1291 cmp_int(fp, "inputrec->eI", -1, ir1->eI, ir2->eI);
1292 cmp_int64(fp, "inputrec->nsteps", ir1->nsteps, ir2->nsteps);
1293 cmp_int64(fp, "inputrec->init_step", ir1->init_step, ir2->init_step);
1294 cmp_int(fp, "inputrec->simulation_part", -1, ir1->simulation_part, ir2->simulation_part);
1295 cmp_int(fp, "inputrec->ePBC", -1, ir1->ePBC, ir2->ePBC);
1296 cmp_bool(fp, "inputrec->bPeriodicMols", -1, ir1->bPeriodicMols, ir2->bPeriodicMols);
1297 cmp_int(fp, "inputrec->cutoff_scheme", -1, ir1->cutoff_scheme, ir2->cutoff_scheme);
1298 cmp_int(fp, "inputrec->nstlist", -1, ir1->nstlist, ir2->nstlist);
1299 cmp_int(fp, "inputrec->nstcomm", -1, ir1->nstcomm, ir2->nstcomm);
1300 cmp_int(fp, "inputrec->comm_mode", -1, ir1->comm_mode, ir2->comm_mode);
1301 cmp_int(fp, "inputrec->nstlog", -1, ir1->nstlog, ir2->nstlog);
1302 cmp_int(fp, "inputrec->nstxout", -1, ir1->nstxout, ir2->nstxout);
1303 cmp_int(fp, "inputrec->nstvout", -1, ir1->nstvout, ir2->nstvout);
1304 cmp_int(fp, "inputrec->nstfout", -1, ir1->nstfout, ir2->nstfout);
1305 cmp_int(fp, "inputrec->nstcalcenergy", -1, ir1->nstcalcenergy, ir2->nstcalcenergy);
1306 cmp_int(fp, "inputrec->nstenergy", -1, ir1->nstenergy, ir2->nstenergy);
1307 cmp_int(fp, "inputrec->nstxout_compressed", -1, ir1->nstxout_compressed, ir2->nstxout_compressed);
1308 cmp_double(fp, "inputrec->init_t", -1, ir1->init_t, ir2->init_t, ftol, abstol);
1309 cmp_double(fp, "inputrec->delta_t", -1, ir1->delta_t, ir2->delta_t, ftol, abstol);
1310 cmp_real(fp, "inputrec->x_compression_precision", -1, ir1->x_compression_precision,
1311 ir2->x_compression_precision, ftol, abstol);
1312 cmp_real(fp, "inputrec->fourierspacing", -1, ir1->fourier_spacing, ir2->fourier_spacing, ftol, abstol);
1313 cmp_int(fp, "inputrec->nkx", -1, ir1->nkx, ir2->nkx);
1314 cmp_int(fp, "inputrec->nky", -1, ir1->nky, ir2->nky);
1315 cmp_int(fp, "inputrec->nkz", -1, ir1->nkz, ir2->nkz);
1316 cmp_int(fp, "inputrec->pme_order", -1, ir1->pme_order, ir2->pme_order);
1317 cmp_real(fp, "inputrec->ewald_rtol", -1, ir1->ewald_rtol, ir2->ewald_rtol, ftol, abstol);
1318 cmp_int(fp, "inputrec->ewald_geometry", -1, ir1->ewald_geometry, ir2->ewald_geometry);
1319 cmp_real(fp, "inputrec->epsilon_surface", -1, ir1->epsilon_surface, ir2->epsilon_surface, ftol, abstol);
1320 cmp_int(fp, "inputrec->bContinuation", -1, static_cast<int>(ir1->bContinuation),
1321 static_cast<int>(ir2->bContinuation));
1322 cmp_int(fp, "inputrec->bShakeSOR", -1, static_cast<int>(ir1->bShakeSOR),
1323 static_cast<int>(ir2->bShakeSOR));
1324 cmp_int(fp, "inputrec->etc", -1, ir1->etc, ir2->etc);
1325 cmp_int(fp, "inputrec->bPrintNHChains", -1, static_cast<int>(ir1->bPrintNHChains),
1326 static_cast<int>(ir2->bPrintNHChains));
1327 cmp_int(fp, "inputrec->epc", -1, ir1->epc, ir2->epc);
1328 cmp_int(fp, "inputrec->epct", -1, ir1->epct, ir2->epct);
1329 cmp_real(fp, "inputrec->tau_p", -1, ir1->tau_p, ir2->tau_p, ftol, abstol);
1330 cmp_rvec(fp, "inputrec->ref_p(x)", -1, ir1->ref_p[XX], ir2->ref_p[XX], ftol, abstol);
1331 cmp_rvec(fp, "inputrec->ref_p(y)", -1, ir1->ref_p[YY], ir2->ref_p[YY], ftol, abstol);
1332 cmp_rvec(fp, "inputrec->ref_p(z)", -1, ir1->ref_p[ZZ], ir2->ref_p[ZZ], ftol, abstol);
1333 cmp_rvec(fp, "inputrec->compress(x)", -1, ir1->compress[XX], ir2->compress[XX], ftol, abstol);
1334 cmp_rvec(fp, "inputrec->compress(y)", -1, ir1->compress[YY], ir2->compress[YY], ftol, abstol);
1335 cmp_rvec(fp, "inputrec->compress(z)", -1, ir1->compress[ZZ], ir2->compress[ZZ], ftol, abstol);
1336 cmp_int(fp, "refcoord_scaling", -1, ir1->refcoord_scaling, ir2->refcoord_scaling);
1337 cmp_rvec(fp, "inputrec->posres_com", -1, ir1->posres_com, ir2->posres_com, ftol, abstol);
1338 cmp_rvec(fp, "inputrec->posres_comB", -1, ir1->posres_comB, ir2->posres_comB, ftol, abstol);
1339 cmp_real(fp, "inputrec->verletbuf_tol", -1, ir1->verletbuf_tol, ir2->verletbuf_tol, ftol, abstol);
1340 cmp_real(fp, "inputrec->rlist", -1, ir1->rlist, ir2->rlist, ftol, abstol);
1341 cmp_real(fp, "inputrec->rtpi", -1, ir1->rtpi, ir2->rtpi, ftol, abstol);
1342 cmp_int(fp, "inputrec->coulombtype", -1, ir1->coulombtype, ir2->coulombtype);
1343 cmp_int(fp, "inputrec->coulomb_modifier", -1, ir1->coulomb_modifier, ir2->coulomb_modifier);
1344 cmp_real(fp, "inputrec->rcoulomb_switch", -1, ir1->rcoulomb_switch, ir2->rcoulomb_switch, ftol, abstol);
1345 cmp_real(fp, "inputrec->rcoulomb", -1, ir1->rcoulomb, ir2->rcoulomb, ftol, abstol);
1346 cmp_int(fp, "inputrec->vdwtype", -1, ir1->vdwtype, ir2->vdwtype);
1347 cmp_int(fp, "inputrec->vdw_modifier", -1, ir1->vdw_modifier, ir2->vdw_modifier);
1348 cmp_real(fp, "inputrec->rvdw_switch", -1, ir1->rvdw_switch, ir2->rvdw_switch, ftol, abstol);
1349 cmp_real(fp, "inputrec->rvdw", -1, ir1->rvdw, ir2->rvdw, ftol, abstol);
1350 cmp_real(fp, "inputrec->epsilon_r", -1, ir1->epsilon_r, ir2->epsilon_r, ftol, abstol);
1351 cmp_real(fp, "inputrec->epsilon_rf", -1, ir1->epsilon_rf, ir2->epsilon_rf, ftol, abstol);
1352 cmp_real(fp, "inputrec->tabext", -1, ir1->tabext, ir2->tabext, ftol, abstol);
1354 cmp_int(fp, "inputrec->eDispCorr", -1, ir1->eDispCorr, ir2->eDispCorr);
1355 cmp_real(fp, "inputrec->shake_tol", -1, ir1->shake_tol, ir2->shake_tol, ftol, abstol);
1356 cmp_int(fp, "inputrec->efep", -1, ir1->efep, ir2->efep);
1357 cmp_fepvals(fp, ir1->fepvals, ir2->fepvals, ftol, abstol);
1358 cmp_int(fp, "inputrec->bSimTemp", -1, static_cast<int>(ir1->bSimTemp), static_cast<int>(ir2->bSimTemp));
1359 if ((ir1->bSimTemp == ir2->bSimTemp) && (ir1->bSimTemp))
1361 cmp_simtempvals(fp, ir1->simtempvals, ir2->simtempvals,
1362 std::min(ir1->fepvals->n_lambda, ir2->fepvals->n_lambda), ftol, abstol);
1364 cmp_int(fp, "inputrec->bExpanded", -1, static_cast<int>(ir1->bExpanded),
1365 static_cast<int>(ir2->bExpanded));
1366 if ((ir1->bExpanded == ir2->bExpanded) && (ir1->bExpanded))
1368 cmp_expandedvals(fp, ir1->expandedvals, ir2->expandedvals,
1369 std::min(ir1->fepvals->n_lambda, ir2->fepvals->n_lambda), ftol, abstol);
1371 cmp_int(fp, "inputrec->nwall", -1, ir1->nwall, ir2->nwall);
1372 cmp_int(fp, "inputrec->wall_type", -1, ir1->wall_type, ir2->wall_type);
1373 cmp_int(fp, "inputrec->wall_atomtype[0]", -1, ir1->wall_atomtype[0], ir2->wall_atomtype[0]);
1374 cmp_int(fp, "inputrec->wall_atomtype[1]", -1, ir1->wall_atomtype[1], ir2->wall_atomtype[1]);
1375 cmp_real(fp, "inputrec->wall_density[0]", -1, ir1->wall_density[0], ir2->wall_density[0], ftol, abstol);
1376 cmp_real(fp, "inputrec->wall_density[1]", -1, ir1->wall_density[1], ir2->wall_density[1], ftol, abstol);
1377 cmp_real(fp, "inputrec->wall_ewald_zfac", -1, ir1->wall_ewald_zfac, ir2->wall_ewald_zfac, ftol, abstol);
1379 cmp_bool(fp, "inputrec->bPull", -1, ir1->bPull, ir2->bPull);
1380 if (ir1->bPull && ir2->bPull)
1385 cmp_bool(fp, "inputrec->bDoAwh", -1, ir1->bDoAwh, ir2->bDoAwh);
1386 if (ir1->bDoAwh && ir2->bDoAwh)
1388 cmp_awhParams(fp, ir1->awhParams, ir2->awhParams, ftol, abstol);
1391 cmp_int(fp, "inputrec->eDisre", -1, ir1->eDisre, ir2->eDisre);
1392 cmp_real(fp, "inputrec->dr_fc", -1, ir1->dr_fc, ir2->dr_fc, ftol, abstol);
1393 cmp_int(fp, "inputrec->eDisreWeighting", -1, ir1->eDisreWeighting, ir2->eDisreWeighting);
1394 cmp_int(fp, "inputrec->bDisreMixed", -1, static_cast<int>(ir1->bDisreMixed),
1395 static_cast<int>(ir2->bDisreMixed));
1396 cmp_int(fp, "inputrec->nstdisreout", -1, ir1->nstdisreout, ir2->nstdisreout);
1397 cmp_real(fp, "inputrec->dr_tau", -1, ir1->dr_tau, ir2->dr_tau, ftol, abstol);
1398 cmp_real(fp, "inputrec->orires_fc", -1, ir1->orires_fc, ir2->orires_fc, ftol, abstol);
1399 cmp_real(fp, "inputrec->orires_tau", -1, ir1->orires_tau, ir2->orires_tau, ftol, abstol);
1400 cmp_int(fp, "inputrec->nstorireout", -1, ir1->nstorireout, ir2->nstorireout);
1401 cmp_real(fp, "inputrec->em_stepsize", -1, ir1->em_stepsize, ir2->em_stepsize, ftol, abstol);
1402 cmp_real(fp, "inputrec->em_tol", -1, ir1->em_tol, ir2->em_tol, ftol, abstol);
1403 cmp_int(fp, "inputrec->niter", -1, ir1->niter, ir2->niter);
1404 cmp_real(fp, "inputrec->fc_stepsize", -1, ir1->fc_stepsize, ir2->fc_stepsize, ftol, abstol);
1405 cmp_int(fp, "inputrec->nstcgsteep", -1, ir1->nstcgsteep, ir2->nstcgsteep);
1406 cmp_int(fp, "inputrec->nbfgscorr", 0, ir1->nbfgscorr, ir2->nbfgscorr);
1407 cmp_int(fp, "inputrec->eConstrAlg", -1, ir1->eConstrAlg, ir2->eConstrAlg);
1408 cmp_int(fp, "inputrec->nProjOrder", -1, ir1->nProjOrder, ir2->nProjOrder);
1409 cmp_real(fp, "inputrec->LincsWarnAngle", -1, ir1->LincsWarnAngle, ir2->LincsWarnAngle, ftol, abstol);
1410 cmp_int(fp, "inputrec->nLincsIter", -1, ir1->nLincsIter, ir2->nLincsIter);
1411 cmp_real(fp, "inputrec->bd_fric", -1, ir1->bd_fric, ir2->bd_fric, ftol, abstol);
1412 cmp_int64(fp, "inputrec->ld_seed", ir1->ld_seed, ir2->ld_seed);
1413 cmp_real(fp, "inputrec->cos_accel", -1, ir1->cos_accel, ir2->cos_accel, ftol, abstol);
1414 cmp_rvec(fp, "inputrec->deform(a)", -1, ir1->deform[XX], ir2->deform[XX], ftol, abstol);
1415 cmp_rvec(fp, "inputrec->deform(b)", -1, ir1->deform[YY], ir2->deform[YY], ftol, abstol);
1416 cmp_rvec(fp, "inputrec->deform(c)", -1, ir1->deform[ZZ], ir2->deform[ZZ], ftol, abstol);
1419 cmp_int(fp, "inputrec->userint1", -1, ir1->userint1, ir2->userint1);
1420 cmp_int(fp, "inputrec->userint2", -1, ir1->userint2, ir2->userint2);
1421 cmp_int(fp, "inputrec->userint3", -1, ir1->userint3, ir2->userint3);
1422 cmp_int(fp, "inputrec->userint4", -1, ir1->userint4, ir2->userint4);
1423 cmp_real(fp, "inputrec->userreal1", -1, ir1->userreal1, ir2->userreal1, ftol, abstol);
1424 cmp_real(fp, "inputrec->userreal2", -1, ir1->userreal2, ir2->userreal2, ftol, abstol);
1425 cmp_real(fp, "inputrec->userreal3", -1, ir1->userreal3, ir2->userreal3, ftol, abstol);
1426 cmp_real(fp, "inputrec->userreal4", -1, ir1->userreal4, ir2->userreal4, ftol, abstol);
1427 cmp_grpopts(fp, &(ir1->opts), &(ir2->opts), ftol, abstol);
1428 gmx::TextWriter writer(fp);
1429 gmx::compareKeyValueTrees(&writer, *ir1->params, *ir2->params, ftol, abstol);
1432 void comp_pull_AB(FILE* fp, pull_params_t* pull, real ftol, real abstol)
1436 for (i = 0; i < pull->ncoord; i++)
1438 fprintf(fp, "comparing pull coord %d\n", i);
1439 cmp_real(fp, "pull-coord->k", -1, pull->coord[i].k, pull->coord[i].kB, ftol, abstol);
1443 gmx_bool inputrecDeform(const t_inputrec* ir)
1445 return (ir->deform[XX][XX] != 0 || ir->deform[YY][YY] != 0 || ir->deform[ZZ][ZZ] != 0
1446 || ir->deform[YY][XX] != 0 || ir->deform[ZZ][XX] != 0 || ir->deform[ZZ][YY] != 0);
1449 gmx_bool inputrecDynamicBox(const t_inputrec* ir)
1451 return (ir->epc != epcNO || ir->eI == eiTPI || inputrecDeform(ir));
1454 gmx_bool inputrecPreserveShape(const t_inputrec* ir)
1456 return (ir->epc != epcNO && ir->deform[XX][XX] == 0
1457 && (ir->epct == epctISOTROPIC || ir->epct == epctSEMIISOTROPIC));
1460 gmx_bool inputrecNeedMutot(const t_inputrec* ir)
1462 return ((ir->coulombtype == eelEWALD || EEL_PME(ir->coulombtype))
1463 && (ir->ewald_geometry == eewg3DC || ir->epsilon_surface != 0));
1466 gmx_bool inputrecExclForces(const t_inputrec* ir)
1468 return (EEL_FULL(ir->coulombtype) || (EEL_RF(ir->coulombtype)));
1471 gmx_bool inputrecNptTrotter(const t_inputrec* ir)
1473 return (((ir->eI == eiVV) || (ir->eI == eiVVAK)) && (ir->epc == epcMTTK) && (ir->etc == etcNOSEHOOVER));
1476 gmx_bool inputrecNvtTrotter(const t_inputrec* ir)
1478 return (((ir->eI == eiVV) || (ir->eI == eiVVAK)) && (ir->epc != epcMTTK) && (ir->etc == etcNOSEHOOVER));
1481 gmx_bool inputrecNphTrotter(const t_inputrec* ir)
1483 return (((ir->eI == eiVV) || (ir->eI == eiVVAK)) && (ir->epc == epcMTTK) && (ir->etc != etcNOSEHOOVER));
1486 bool inputrecPbcXY2Walls(const t_inputrec* ir)
1488 return (ir->ePBC == epbcXY && ir->nwall == 2);
1491 bool integratorHasConservedEnergyQuantity(const t_inputrec* ir)
1495 // Energy minimization or stochastic integrator: no conservation
1498 else if (ir->etc == etcNO && ir->epc == epcNO)
1500 // The total energy is conserved, no additional conserved quanitity
1505 // Shear stress with Parrinello-Rahman is not supported (tedious)
1507 ((ir->epc == epcPARRINELLORAHMAN || ir->epc == epcMTTK)
1508 && (ir->ref_p[YY][XX] != 0 || ir->ref_p[ZZ][XX] != 0 || ir->ref_p[ZZ][YY] != 0));
1510 return !ETC_ANDERSEN(ir->etc) && !shearWithPR;
1514 bool integratorHasReferenceTemperature(const t_inputrec* ir)
1516 return ((ir->etc != etcNO) || EI_SD(ir->eI) || (ir->eI == eiBD) || EI_TPI(ir->eI));
1519 int inputrec2nboundeddim(const t_inputrec* ir)
1521 if (inputrecPbcXY2Walls(ir))
1527 return ePBC2npbcdim(ir->ePBC);
1531 int ndof_com(const t_inputrec* ir)
1538 case epbcNONE: n = 3; break;
1539 case epbcXY: n = (ir->nwall == 0 ? 3 : 2); break;
1540 case epbcSCREW: n = 1; break;
1541 default: gmx_incons("Unknown pbc in calc_nrdf");
1547 real maxReferenceTemperature(const t_inputrec& ir)
1549 if (EI_ENERGY_MINIMIZATION(ir.eI) || ir.eI == eiNM)
1554 if (EI_MD(ir.eI) && ir.etc == etcNO)
1559 /* SD and BD also use ref_t and tau_t for setting the reference temperature.
1560 * TPI can be treated as MD, since it needs an ensemble temperature.
1563 real maxTemperature = 0;
1564 for (int i = 0; i < ir.opts.ngtc; i++)
1566 if (ir.opts.tau_t[i] >= 0)
1568 maxTemperature = std::max(maxTemperature, ir.opts.ref_t[i]);
1572 return maxTemperature;
1575 bool haveEwaldSurfaceContribution(const t_inputrec& ir)
1577 return EEL_PME_EWALD(ir.coulombtype) && (ir.ewald_geometry == eewg3DC || ir.epsilon_surface != 0);