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 by the GROMACS development team.
7 * Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by
8 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
9 * and including many others, as listed in the AUTHORS file in the
10 * top-level source directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
40 #include "gromacs/utility/enumerationhelpers.h"
51 #include "gromacs/applied_forces/awh/read_params.h"
52 #include "gromacs/math/veccompare.h"
53 #include "gromacs/math/vecdump.h"
54 #include "gromacs/mdlib/vcm.h"
55 #include "gromacs/mdtypes/awh_params.h"
56 #include "gromacs/mdtypes/md_enums.h"
57 #include "gromacs/mdtypes/multipletimestepping.h"
58 #include "gromacs/mdtypes/pull_params.h"
59 #include "gromacs/pbcutil/pbc.h"
60 #include "gromacs/utility/compare.h"
61 #include "gromacs/utility/cstringutil.h"
62 #include "gromacs/utility/fatalerror.h"
63 #include "gromacs/utility/keyvaluetree.h"
64 #include "gromacs/utility/smalloc.h"
65 #include "gromacs/utility/snprintf.h"
66 #include "gromacs/utility/strconvert.h"
67 #include "gromacs/utility/stringutil.h"
68 #include "gromacs/utility/textwriter.h"
69 #include "gromacs/utility/txtdump.h"
71 //! Macro to select a bool name
72 #define EBOOL(e) gmx::boolToString(e)
74 /* Default values for nstcalcenergy, used when the are no other restrictions. */
75 constexpr int c_defaultNstCalcEnergy = 10;
77 /* The minimum number of integration steps required for reasonably accurate
78 * integration of first and second order coupling algorithms.
80 const int nstmin_berendsen_tcouple = 5;
81 const int nstmin_berendsen_pcouple = 10;
82 const int nstmin_harmonic = 20;
84 /* Default values for T- and P- coupling intervals, used when the are no other
87 constexpr int c_defaultNstTCouple = 10;
88 constexpr int c_defaultNstPCouple = 10;
90 t_inputrec::t_inputrec() :
91 fepvals(std::make_unique<t_lambda>()),
92 simtempvals(std::make_unique<t_simtemp>()),
93 expandedvals(std::make_unique<t_expanded>())
97 t_inputrec::~t_inputrec()
102 int ir_optimal_nstcalcenergy(const t_inputrec* ir)
112 nst = c_defaultNstCalcEnergy;
117 nst = std::lcm(nst, ir->mtsLevels.back().stepFactor);
123 int tcouple_min_integration_steps(TemperatureCoupling etc)
129 case TemperatureCoupling::No: n = 0; break;
130 case TemperatureCoupling::Berendsen:
131 case TemperatureCoupling::Yes: n = nstmin_berendsen_tcouple; break;
132 case TemperatureCoupling::VRescale:
133 /* V-rescale supports instantaneous rescaling */
136 case TemperatureCoupling::NoseHoover: n = nstmin_harmonic; break;
137 case TemperatureCoupling::Andersen:
138 case TemperatureCoupling::AndersenMassive: n = 1; break;
139 default: gmx_incons("Unknown etc value");
145 int ir_optimal_nsttcouple(const t_inputrec* ir)
147 int nmin, nwanted, n;
151 nmin = tcouple_min_integration_steps(ir->etc);
153 nwanted = c_defaultNstTCouple;
156 if (ir->etc != TemperatureCoupling::No)
158 for (g = 0; g < ir->opts.ngtc; g++)
160 if (ir->opts.tau_t[g] > 0)
162 tau_min = std::min(tau_min, ir->opts.tau_t[g]);
167 if (nmin == 0 || ir->delta_t * nwanted <= tau_min)
173 n = static_cast<int>(tau_min / (ir->delta_t * nmin) + 0.001);
178 while (nwanted % n != 0)
187 int pcouple_min_integration_steps(PressureCoupling epc)
193 case PressureCoupling::No: n = 0; break;
194 case PressureCoupling::Berendsen:
195 case PressureCoupling::CRescale:
196 case PressureCoupling::Isotropic: n = nstmin_berendsen_pcouple; break;
197 case PressureCoupling::ParrinelloRahman:
198 case PressureCoupling::Mttk: n = nstmin_harmonic; break;
199 default: gmx_incons("Unknown epc value");
205 int ir_optimal_nstpcouple(const t_inputrec* ir)
207 const int minIntegrationSteps = pcouple_min_integration_steps(ir->epc);
209 const int nwanted = c_defaultNstPCouple;
211 // With multiple time-stepping we can only compute the pressure at slowest steps
212 const int minNstPCouple = (ir->useMts ? ir->mtsLevels.back().stepFactor : 1);
215 if (minIntegrationSteps == 0 || ir->delta_t * nwanted <= ir->tau_p)
221 n = static_cast<int>(ir->tau_p / (ir->delta_t * minIntegrationSteps) + 0.001);
222 if (n < minNstPCouple)
226 // Without MTS we try to make nstpcouple a "nice" number
229 while (nwanted % n != 0)
236 // With MTS, nstpcouple should be a multiple of the slowest MTS interval
239 n = n - (n % minNstPCouple);
245 gmx_bool ir_coulomb_switched(const t_inputrec* ir)
247 return (ir->coulombtype == CoulombInteractionType::Switch
248 || ir->coulombtype == CoulombInteractionType::Shift
249 || ir->coulombtype == CoulombInteractionType::PmeSwitch
250 || ir->coulombtype == CoulombInteractionType::PmeUserSwitch
251 || ir->coulomb_modifier == InteractionModifiers::PotSwitch
252 || ir->coulomb_modifier == InteractionModifiers::ForceSwitch);
255 gmx_bool ir_coulomb_is_zero_at_cutoff(const t_inputrec* ir)
257 return (ir->cutoff_scheme == CutoffScheme::Verlet || ir_coulomb_switched(ir)
258 || ir->coulomb_modifier != InteractionModifiers::None
259 || ir->coulombtype == CoulombInteractionType::RFZero);
262 gmx_bool ir_coulomb_might_be_zero_at_cutoff(const t_inputrec* ir)
264 return (ir_coulomb_is_zero_at_cutoff(ir) || ir->coulombtype == CoulombInteractionType::User
265 || ir->coulombtype == CoulombInteractionType::PmeUser);
268 gmx_bool ir_vdw_switched(const t_inputrec* ir)
270 return (ir->vdwtype == VanDerWaalsType::Switch || ir->vdwtype == VanDerWaalsType::Shift
271 || ir->vdw_modifier == InteractionModifiers::PotSwitch
272 || ir->vdw_modifier == InteractionModifiers::ForceSwitch);
275 gmx_bool ir_vdw_is_zero_at_cutoff(const t_inputrec* ir)
277 return (ir->cutoff_scheme == CutoffScheme::Verlet || ir_vdw_switched(ir)
278 || ir->vdw_modifier != InteractionModifiers::None);
281 gmx_bool ir_vdw_might_be_zero_at_cutoff(const t_inputrec* ir)
283 return (ir_vdw_is_zero_at_cutoff(ir) || ir->vdwtype == VanDerWaalsType::User);
286 static void done_t_rot(t_rot* rot)
292 if (rot->grp != nullptr)
294 for (int i = 0; i < rot->ngrp; i++)
296 sfree(rot->grp[i].ind);
297 sfree(rot->grp[i].x_ref);
304 static void done_t_swapCoords(t_swapcoords* swapCoords)
306 if (swapCoords == nullptr)
310 for (int i = 0; i < swapCoords->ngrp; i++)
312 sfree(swapCoords->grp[i].ind);
313 sfree(swapCoords->grp[i].molname);
315 sfree(swapCoords->grp);
319 void done_inputrec(t_inputrec* ir)
321 sfree(ir->opts.nrdf);
322 sfree(ir->opts.ref_t);
323 for (int i = 0; i < ir->opts.ngtc; i++)
325 sfree(ir->opts.anneal_time[i]);
326 sfree(ir->opts.anneal_temp[i]);
328 sfree(ir->opts.annealing);
329 sfree(ir->opts.anneal_npoints);
330 sfree(ir->opts.anneal_time);
331 sfree(ir->opts.anneal_temp);
332 sfree(ir->opts.tau_t);
333 sfree(ir->opts.acceleration);
334 sfree(ir->opts.nFreeze);
335 sfree(ir->opts.egp_flags);
337 done_t_swapCoords(ir->swap);
342 static void pr_grp_opts(FILE* out, int indent, const char* title, const t_grpopts* opts, gmx_bool bMDPformat)
348 fprintf(out, "%s:\n", title);
351 pr_indent(out, indent);
352 fprintf(out, "nrdf%s", bMDPformat ? " = " : ":");
353 for (i = 0; (i < opts->ngtc); i++)
355 fprintf(out, " %10g", opts->nrdf[i]);
359 pr_indent(out, indent);
360 fprintf(out, "ref-t%s", bMDPformat ? " = " : ":");
361 for (i = 0; (i < opts->ngtc); i++)
363 fprintf(out, " %10g", opts->ref_t[i]);
367 pr_indent(out, indent);
368 fprintf(out, "tau-t%s", bMDPformat ? " = " : ":");
369 for (i = 0; (i < opts->ngtc); i++)
371 fprintf(out, " %10g", opts->tau_t[i]);
375 /* Pretty-print the simulated annealing info */
376 fprintf(out, "annealing%s", bMDPformat ? " = " : ":");
377 for (i = 0; (i < opts->ngtc); i++)
379 fprintf(out, " %10s", enumValueToString(opts->annealing[i]));
383 fprintf(out, "annealing-npoints%s", bMDPformat ? " = " : ":");
384 for (i = 0; (i < opts->ngtc); i++)
386 fprintf(out, " %10d", opts->anneal_npoints[i]);
390 for (i = 0; (i < opts->ngtc); i++)
392 if (opts->anneal_npoints[i] > 0)
394 fprintf(out, "annealing-time [%d]:\t", i);
395 for (j = 0; (j < opts->anneal_npoints[i]); j++)
397 fprintf(out, " %10.1f", opts->anneal_time[i][j]);
400 fprintf(out, "annealing-temp [%d]:\t", i);
401 for (j = 0; (j < opts->anneal_npoints[i]); j++)
403 fprintf(out, " %10.1f", opts->anneal_temp[i][j]);
409 pr_indent(out, indent);
410 fprintf(out, "acc:\t");
411 for (i = 0; (i < opts->ngacc); i++)
413 for (m = 0; (m < DIM); m++)
415 fprintf(out, " %10g", opts->acceleration[i][m]);
420 pr_indent(out, indent);
421 fprintf(out, "nfreeze:");
422 for (i = 0; (i < opts->ngfrz); i++)
424 for (m = 0; (m < DIM); m++)
426 fprintf(out, " %10s", opts->nFreeze[i][m] ? "Y" : "N");
432 for (i = 0; (i < opts->ngener); i++)
434 pr_indent(out, indent);
435 fprintf(out, "energygrp-flags[%3d]:", i);
436 for (m = 0; (m < opts->ngener); m++)
438 fprintf(out, " %d", opts->egp_flags[opts->ngener * i + m]);
446 static void pr_matrix(FILE* fp, int indent, const char* title, const rvec* m, gmx_bool bMDPformat)
451 "%-10s = %g %g %g %g %g %g\n",
462 pr_rvecs(fp, indent, title, m, DIM);
466 #define PS(t, s) pr_str(fp, indent, t, s)
467 #define PI(t, s) pr_int(fp, indent, t, s)
468 #define PSTEP(t, s) pr_int64(fp, indent, t, s)
469 #define PR(t, s) pr_real(fp, indent, t, s)
470 #define PD(t, s) pr_double(fp, indent, t, s)
472 static void pr_pull_group(FILE* fp, int indent, int g, const t_pull_group* pgrp)
474 pr_indent(fp, indent);
475 fprintf(fp, "pull-group %d:\n", g);
477 pr_ivec_block(fp, indent, "atom", pgrp->ind.data(), pgrp->ind.size(), TRUE);
478 pr_rvec(fp, indent, "weight", pgrp->weight.data(), pgrp->weight.size(), TRUE);
479 PI("pbcatom", pgrp->pbcatom);
482 static void pr_pull_coord(FILE* fp, int indent, int c, const t_pull_coord* pcrd)
484 pr_indent(fp, indent);
485 fprintf(fp, "pull-coord %d:\n", c);
486 PS("type", enumValueToString(pcrd->eType));
487 if (pcrd->eType == PullingAlgorithm::External)
489 PS("potential-provider", pcrd->externalPotentialProvider.c_str());
491 PS("geometry", enumValueToString(pcrd->eGeom));
492 for (int g = 0; g < pcrd->ngroup; g++)
494 std::string buffer = gmx::formatString("group[%d]", g);
495 PI(buffer.c_str(), pcrd->group[g]);
497 pr_ivec(fp, indent, "dim", pcrd->dim, DIM, TRUE);
498 pr_rvec(fp, indent, "origin", pcrd->origin, DIM, TRUE);
499 pr_rvec(fp, indent, "vec", pcrd->vec, DIM, TRUE);
500 PS("start", EBOOL(pcrd->bStart));
501 PR("init", pcrd->init);
502 PR("rate", pcrd->rate);
507 static void pr_simtempvals(FILE* fp, int indent, const t_simtemp* simtemp, int n_lambda)
509 PS("simulated-tempering-scaling", enumValueToString(simtemp->eSimTempScale));
510 PR("sim-temp-low", simtemp->simtemp_low);
511 PR("sim-temp-high", simtemp->simtemp_high);
512 pr_rvec(fp, indent, "simulated tempering temperatures", simtemp->temperatures.data(), n_lambda, TRUE);
515 static void pr_expandedvals(FILE* fp, int indent, const t_expanded* expand, int n_lambda)
518 PI("nstexpanded", expand->nstexpanded);
519 PS("lmc-stats", enumValueToString(expand->elamstats));
520 PS("lmc-move", enumValueToString(expand->elmcmove));
521 PS("lmc-weights-equil", enumValueToString(expand->elmceq));
522 if (expand->elmceq == LambdaWeightWillReachEquilibrium::NumAtLambda)
524 PI("weight-equil-number-all-lambda", expand->equil_n_at_lam);
526 if (expand->elmceq == LambdaWeightWillReachEquilibrium::Samples)
528 PI("weight-equil-number-samples", expand->equil_samples);
530 if (expand->elmceq == LambdaWeightWillReachEquilibrium::Steps)
532 PI("weight-equil-number-steps", expand->equil_steps);
534 if (expand->elmceq == LambdaWeightWillReachEquilibrium::WLDelta)
536 PR("weight-equil-wl-delta", expand->equil_wl_delta);
538 if (expand->elmceq == LambdaWeightWillReachEquilibrium::Ratio)
540 PR("weight-equil-count-ratio", expand->equil_ratio);
542 PI("lmc-seed", expand->lmc_seed);
543 PR("mc-temperature", expand->mc_temp);
544 PI("lmc-repeats", expand->lmc_repeats);
545 PI("lmc-gibbsdelta", expand->gibbsdeltalam);
546 PI("lmc-forced-nstart", expand->lmc_forced_nstart);
547 PS("symmetrized-transition-matrix", EBOOL(expand->bSymmetrizedTMatrix));
548 PI("nst-transition-matrix", expand->nstTij);
549 PI("mininum-var-min", expand->minvarmin); /*default is reasonable */
550 PI("weight-c-range", expand->c_range); /* default is just C=0 */
551 PR("wl-scale", expand->wl_scale);
552 PR("wl-ratio", expand->wl_ratio);
553 PR("init-wl-delta", expand->init_wl_delta);
554 PS("wl-oneovert", EBOOL(expand->bWLoneovert));
556 pr_indent(fp, indent);
557 pr_rvec(fp, indent, "init-lambda-weights", expand->init_lambda_weights.data(), n_lambda, TRUE);
558 PS("init-weights", EBOOL(expand->bInit_weights));
561 static void pr_fepvals(FILE* fp, int indent, const t_lambda* fep, gmx_bool bMDPformat)
565 PR("init-lambda", fep->init_lambda);
566 PI("init-lambda-state", fep->init_fep_state);
567 PR("delta-lambda", fep->delta_lambda);
568 PI("nstdhdl", fep->nstdhdl);
572 PI("n-lambdas", fep->n_lambda);
574 if (fep->n_lambda > 0)
576 pr_indent(fp, indent);
577 fprintf(fp, "separate-dvdl%s\n", bMDPformat ? " = " : ":");
578 for (auto i : gmx::EnumerationArray<FreeEnergyPerturbationCouplingType, bool>::keys())
580 fprintf(fp, "%18s = ", enumValueToString(i));
581 if (fep->separate_dvdl[i])
583 fprintf(fp, " TRUE");
587 fprintf(fp, " FALSE");
591 fprintf(fp, "all-lambdas%s\n", bMDPformat ? " = " : ":");
592 for (auto key : gmx::EnumerationArray<FreeEnergyPerturbationCouplingType, bool>::keys())
594 fprintf(fp, "%18s = ", enumValueToString(key));
595 int i = static_cast<int>(key);
596 for (j = 0; j < fep->n_lambda; j++)
598 fprintf(fp, " %10g", fep->all_lambda[i][j]);
603 PI("calc-lambda-neighbors", fep->lambda_neighbors);
604 PS("dhdl-print-energy", enumValueToString(fep->edHdLPrintEnergy));
605 PR("sc-alpha", fep->sc_alpha);
606 PI("sc-power", fep->sc_power);
607 PR("sc-r-power", fep->sc_r_power);
608 PR("sc-sigma", fep->sc_sigma);
609 PR("sc-sigma-min", fep->sc_sigma_min);
610 PS("sc-coul", EBOOL(fep->bScCoul));
611 PI("dh-hist-size", fep->dh_hist_size);
612 PD("dh-hist-spacing", fep->dh_hist_spacing);
613 PS("separate-dhdl-file", enumValueToString(fep->separate_dhdl_file));
614 PS("dhdl-derivatives", enumValueToString(fep->dhdl_derivatives));
615 PS("sc-function", enumValueToString(fep->softcoreFunction));
616 PR("sc-gapsys-scale-linpoint-lj", fep->scGapsysScaleLinpointLJ);
617 PR("sc-gapsys-scale-linpoint-q", fep->scGapsysScaleLinpointQ);
618 PR("sc-gapsys-sigma-lj", fep->scGapsysSigmaLJ);
621 static void pr_pull(FILE* fp, int indent, const pull_params_t& pull)
625 PR("pull-cylinder-r", pull.cylinder_r);
626 PR("pull-constr-tol", pull.constr_tol);
627 PS("pull-print-COM", EBOOL(pull.bPrintCOM));
628 PS("pull-print-ref-value", EBOOL(pull.bPrintRefValue));
629 PS("pull-print-components", EBOOL(pull.bPrintComp));
630 PI("pull-nstxout", pull.nstxout);
631 PI("pull-nstfout", pull.nstfout);
632 PS("pull-pbc-ref-prev-step-com", EBOOL(pull.bSetPbcRefToPrevStepCOM));
633 PS("pull-xout-average", EBOOL(pull.bXOutAverage));
634 PS("pull-fout-average", EBOOL(pull.bFOutAverage));
635 PI("pull-ngroups", pull.ngroup);
636 for (g = 0; g < pull.ngroup; g++)
638 pr_pull_group(fp, indent, g, &pull.group[g]);
640 PI("pull-ncoords", pull.ncoord);
641 for (g = 0; g < pull.ncoord; g++)
643 pr_pull_coord(fp, indent, g, &pull.coord[g]);
647 static void pr_awh_bias_dim(FILE* fp, int indent, const gmx::AwhDimParams& awhDimParams, const char* prefix)
649 pr_indent(fp, indent);
651 fprintf(fp, "%s:\n", prefix);
652 PS("coord-provider", enumValueToString(awhDimParams.coordinateProvider()));
653 PI("coord-index", awhDimParams.coordinateIndex() + 1);
654 PR("start", awhDimParams.origin());
655 PR("end", awhDimParams.end());
656 PR("period", awhDimParams.period());
657 PR("force-constant", awhDimParams.forceConstant());
658 PR("diffusion", awhDimParams.diffusion());
659 PR("cover-diameter", awhDimParams.coverDiameter());
662 static void pr_awh_bias(FILE* fp, int indent, const gmx::AwhBiasParams& awhBiasParams, const char* prefix)
666 sprintf(opt, "%s-error-init", prefix);
667 PR(opt, awhBiasParams.initialErrorEstimate());
668 sprintf(opt, "%s-growth", prefix);
669 PS(opt, enumValueToString(awhBiasParams.growthType()));
670 sprintf(opt, "%s-target", prefix);
671 PS(opt, enumValueToString(awhBiasParams.targetDistribution()));
672 sprintf(opt, "%s-target-beta-scalng", prefix);
673 PR(opt, awhBiasParams.targetBetaScaling());
674 sprintf(opt, "%s-target-cutoff", prefix);
675 PR(opt, awhBiasParams.targetCutoff());
676 sprintf(opt, "%s-user-data", prefix);
677 PS(opt, EBOOL(awhBiasParams.userPMFEstimate()));
678 sprintf(opt, "%s-share-group", prefix);
679 PI(opt, awhBiasParams.shareGroup());
680 sprintf(opt, "%s-equilibrate-histogram", prefix);
681 PS(opt, EBOOL(awhBiasParams.equilibrateHistogram()));
682 sprintf(opt, "%s-ndim", prefix);
683 PI(opt, awhBiasParams.ndim());
686 for (const auto& dimParam : awhBiasParams.dimParams())
688 char prefixdim[STRLEN];
689 sprintf(prefixdim, "%s-dim%d", prefix, d + 1);
690 pr_awh_bias_dim(fp, indent, dimParam, prefixdim);
695 static void pr_awh(FILE* fp, int indent, gmx::AwhParams* awhParams)
697 PS("awh-potential", enumValueToString(awhParams->potential()));
698 PI("awh-seed", awhParams->seed());
699 PI("awh-nstout", awhParams->nstout());
700 PI("awh-nstsample", awhParams->nstSampleCoord());
701 PI("awh-nsamples-update", awhParams->numSamplesUpdateFreeEnergy());
702 PS("awh-share-bias-multisim", EBOOL(awhParams->shareBiasMultisim()));
703 PI("awh-nbias", awhParams->numBias());
706 for (const auto& awhBiasParam : awhParams->awhBiasParams())
708 auto prefix = gmx::formatString("awh%d", k + 1);
709 pr_awh_bias(fp, indent, awhBiasParam, prefix.c_str());
714 static void pr_rotgrp(FILE* fp, int indent, int g, const t_rotgrp* rotg)
716 pr_indent(fp, indent);
717 fprintf(fp, "rot-group %d:\n", g);
719 PS("rot-type", enumValueToString(rotg->eType));
720 PS("rot-massw", EBOOL(rotg->bMassW));
721 pr_ivec_block(fp, indent, "atom", rotg->ind, rotg->nat, TRUE);
722 pr_rvecs(fp, indent, "x-ref", rotg->x_ref, rotg->nat);
723 pr_rvec(fp, indent, "rot-vec", rotg->inputVec, DIM, TRUE);
724 pr_rvec(fp, indent, "rot-pivot", rotg->pivot, DIM, TRUE);
725 PR("rot-rate", rotg->rate);
726 PR("rot-k", rotg->k);
727 PR("rot-slab-dist", rotg->slab_dist);
728 PR("rot-min-gauss", rotg->min_gaussian);
729 PR("rot-eps", rotg->eps);
730 PS("rot-fit-method", enumValueToString(rotg->eFittype));
731 PI("rot-potfit-nstep", rotg->PotAngle_nstep);
732 PR("rot-potfit-step", rotg->PotAngle_step);
735 static void pr_rot(FILE* fp, int indent, const t_rot* rot)
739 PI("rot-nstrout", rot->nstrout);
740 PI("rot-nstsout", rot->nstsout);
741 PI("rot-ngroups", rot->ngrp);
742 for (g = 0; g < rot->ngrp; g++)
744 pr_rotgrp(fp, indent, g, &rot->grp[g]);
749 static void pr_swap(FILE* fp, int indent, const t_swapcoords* swap)
753 /* Enums for better readability of the code */
761 PI("swap-frequency", swap->nstswap);
763 /* The split groups that define the compartments */
764 for (int j = 0; j < 2; j++)
766 snprintf(str, STRLEN, "massw_split%d", j);
767 PS(str, EBOOL(swap->massw_split[j]));
768 snprintf(str, STRLEN, "split atoms group %d", j);
769 pr_ivec_block(fp, indent, str, swap->grp[j].ind, swap->grp[j].nat, TRUE);
772 /* The solvent group */
776 swap->grp[static_cast<int>(SwapGroupSplittingType::Solvent)].molname);
780 swap->grp[static_cast<int>(SwapGroupSplittingType::Solvent)].ind,
781 swap->grp[static_cast<int>(SwapGroupSplittingType::Solvent)].nat,
784 /* Now print the indices for all the ion groups: */
785 for (int ig = static_cast<int>(SwapGroupSplittingType::Count); ig < swap->ngrp; ig++)
787 snprintf(str, STRLEN, "ion group %s", swap->grp[ig].molname);
788 pr_ivec_block(fp, indent, str, swap->grp[ig].ind, swap->grp[ig].nat, TRUE);
791 PR("cyl0-r", swap->cyl0r);
792 PR("cyl0-up", swap->cyl0u);
793 PR("cyl0-down", swap->cyl0l);
794 PR("cyl1-r", swap->cyl1r);
795 PR("cyl1-up", swap->cyl1u);
796 PR("cyl1-down", swap->cyl1l);
797 PI("coupl-steps", swap->nAverage);
799 /* Print the requested ion counts for both compartments */
800 for (int ic = eCompA; ic <= eCompB; ic++)
802 for (int ig = static_cast<int>(SwapGroupSplittingType::Count); ig < swap->ngrp; ig++)
804 snprintf(str, STRLEN, "%s-in-%c", swap->grp[ig].molname, 'A' + ic);
805 PI(str, swap->grp[ig].nmolReq[ic]);
809 PR("threshold", swap->threshold);
810 PR("bulk-offsetA", swap->bulkOffset[eCompA]);
811 PR("bulk-offsetB", swap->bulkOffset[eCompB]);
815 static void pr_imd(FILE* fp, int indent, const t_IMD* imd)
817 PI("IMD-atoms", imd->nat);
818 pr_ivec_block(fp, indent, "atom", imd->ind, imd->nat, TRUE);
822 void pr_inputrec(FILE* fp, int indent, const char* title, const t_inputrec* ir, gmx_bool bMDPformat)
824 const char* infbuf = "inf";
826 if (available(fp, ir, indent, title))
830 indent = pr_title(fp, indent, title);
832 /* Try to make this list appear in the same order as the
833 * options are written in the default mdout.mdp, and with
834 * the same user-exposed names to facilitate debugging.
836 PS("integrator", enumValueToString(ir->eI));
837 PR("tinit", ir->init_t);
838 PR("dt", ir->delta_t);
839 PSTEP("nsteps", ir->nsteps);
840 PSTEP("init-step", ir->init_step);
841 PI("simulation-part", ir->simulation_part);
842 PS("mts", EBOOL(ir->useMts));
845 for (int mtsIndex = 1; mtsIndex < static_cast<int>(ir->mtsLevels.size()); mtsIndex++)
847 const auto& mtsLevel = ir->mtsLevels[mtsIndex];
848 const std::string forceKey = gmx::formatString("mts-level%d-forces", mtsIndex + 1);
849 std::string forceGroups;
850 for (int i = 0; i < static_cast<int>(gmx::MtsForceGroups::Count); i++)
852 if (mtsLevel.forceGroups[i])
854 if (!forceGroups.empty())
858 forceGroups += gmx::mtsForceGroupNames[i];
861 PS(forceKey.c_str(), forceGroups.c_str());
862 const std::string factorKey = gmx::formatString("mts-level%d-factor", mtsIndex + 1);
863 PI(factorKey.c_str(), mtsLevel.stepFactor);
866 PS("comm-mode", enumValueToString(ir->comm_mode));
867 PI("nstcomm", ir->nstcomm);
869 /* Langevin dynamics */
870 PR("bd-fric", ir->bd_fric);
871 PSTEP("ld-seed", ir->ld_seed);
873 /* Energy minimization */
874 PR("emtol", ir->em_tol);
875 PR("emstep", ir->em_stepsize);
876 PI("niter", ir->niter);
877 PR("fcstep", ir->fc_stepsize);
878 PI("nstcgsteep", ir->nstcgsteep);
879 PI("nbfgscorr", ir->nbfgscorr);
881 /* Test particle insertion */
882 PR("rtpi", ir->rtpi);
885 PI("nstxout", ir->nstxout);
886 PI("nstvout", ir->nstvout);
887 PI("nstfout", ir->nstfout);
888 PI("nstlog", ir->nstlog);
889 PI("nstcalcenergy", ir->nstcalcenergy);
890 PI("nstenergy", ir->nstenergy);
891 PI("nstxout-compressed", ir->nstxout_compressed);
892 PR("compressed-x-precision", ir->x_compression_precision);
894 /* Neighborsearching parameters */
895 PS("cutoff-scheme", enumValueToString(ir->cutoff_scheme));
896 PI("nstlist", ir->nstlist);
897 PS("pbc", c_pbcTypeNames[ir->pbcType].c_str());
898 PS("periodic-molecules", EBOOL(ir->bPeriodicMols));
899 PR("verlet-buffer-tolerance", ir->verletbuf_tol);
900 PR("rlist", ir->rlist);
902 /* Options for electrostatics and VdW */
903 PS("coulombtype", enumValueToString(ir->coulombtype));
904 PS("coulomb-modifier", enumValueToString(ir->coulomb_modifier));
905 PR("rcoulomb-switch", ir->rcoulomb_switch);
906 PR("rcoulomb", ir->rcoulomb);
907 if (ir->epsilon_r != 0)
909 PR("epsilon-r", ir->epsilon_r);
913 PS("epsilon-r", infbuf);
915 if (ir->epsilon_rf != 0)
917 PR("epsilon-rf", ir->epsilon_rf);
921 PS("epsilon-rf", infbuf);
923 PS("vdw-type", enumValueToString(ir->vdwtype));
924 PS("vdw-modifier", enumValueToString(ir->vdw_modifier));
925 PR("rvdw-switch", ir->rvdw_switch);
926 PR("rvdw", ir->rvdw);
927 PS("DispCorr", enumValueToString(ir->eDispCorr));
928 PR("table-extension", ir->tabext);
930 PR("fourierspacing", ir->fourier_spacing);
931 PI("fourier-nx", ir->nkx);
932 PI("fourier-ny", ir->nky);
933 PI("fourier-nz", ir->nkz);
934 PI("pme-order", ir->pme_order);
935 PR("ewald-rtol", ir->ewald_rtol);
936 PR("ewald-rtol-lj", ir->ewald_rtol_lj);
937 PS("lj-pme-comb-rule", enumValueToString(ir->ljpme_combination_rule));
938 PS("ewald-geometry", enumValueToString(ir->ewald_geometry));
939 PR("epsilon-surface", ir->epsilon_surface);
941 /* Options for weak coupling algorithms */
942 PS("tcoupl", enumValueToString(ir->etc));
943 PI("nsttcouple", ir->nsttcouple);
944 PI("nh-chain-length", ir->opts.nhchainlength);
945 PS("print-nose-hoover-chain-variables", EBOOL(ir->bPrintNHChains));
947 PS("pcoupl", enumValueToString(ir->epc));
948 PS("pcoupltype", enumValueToString(ir->epct));
949 PI("nstpcouple", ir->nstpcouple);
950 PR("tau-p", ir->tau_p);
951 pr_matrix(fp, indent, "compressibility", ir->compress, bMDPformat);
952 pr_matrix(fp, indent, "ref-p", ir->ref_p, bMDPformat);
953 PS("refcoord-scaling", enumValueToString(ir->refcoord_scaling));
957 fprintf(fp, "posres-com = %g %g %g\n", ir->posres_com[XX], ir->posres_com[YY], ir->posres_com[ZZ]);
958 fprintf(fp, "posres-comB = %g %g %g\n", ir->posres_comB[XX], ir->posres_comB[YY], ir->posres_comB[ZZ]);
962 pr_rvec(fp, indent, "posres-com", ir->posres_com, DIM, TRUE);
963 pr_rvec(fp, indent, "posres-comB", ir->posres_comB, DIM, TRUE);
967 PS("QMMM", EBOOL(ir->bQMMM));
968 fprintf(fp, "%s:\n", "qm-opts");
969 pr_int(fp, indent, "ngQM", ir->opts.ngQM);
971 /* CONSTRAINT OPTIONS */
972 PS("constraint-algorithm", enumValueToString(ir->eConstrAlg));
973 PS("continuation", EBOOL(ir->bContinuation));
975 PS("Shake-SOR", EBOOL(ir->bShakeSOR));
976 PR("shake-tol", ir->shake_tol);
977 PI("lincs-order", ir->nProjOrder);
978 PI("lincs-iter", ir->nLincsIter);
979 PR("lincs-warnangle", ir->LincsWarnAngle);
982 PI("nwall", ir->nwall);
983 PS("wall-type", enumValueToString(ir->wall_type));
984 PR("wall-r-linpot", ir->wall_r_linpot);
986 PI("wall-atomtype[0]", ir->wall_atomtype[0]);
987 PI("wall-atomtype[1]", ir->wall_atomtype[1]);
989 PR("wall-density[0]", ir->wall_density[0]);
990 PR("wall-density[1]", ir->wall_density[1]);
991 PR("wall-ewald-zfac", ir->wall_ewald_zfac);
994 PS("pull", EBOOL(ir->bPull));
997 pr_pull(fp, indent, *ir->pull);
1001 PS("awh", EBOOL(ir->bDoAwh));
1004 pr_awh(fp, indent, ir->awhParams.get());
1007 /* ENFORCED ROTATION */
1008 PS("rotation", EBOOL(ir->bRot));
1011 pr_rot(fp, indent, ir->rot);
1014 /* INTERACTIVE MD */
1015 PS("interactiveMD", EBOOL(ir->bIMD));
1018 pr_imd(fp, indent, ir->imd);
1021 /* NMR refinement stuff */
1022 PS("disre", enumValueToString(ir->eDisre));
1023 PS("disre-weighting", enumValueToString(ir->eDisreWeighting));
1024 PS("disre-mixed", EBOOL(ir->bDisreMixed));
1025 PR("dr-fc", ir->dr_fc);
1026 PR("dr-tau", ir->dr_tau);
1027 PR("nstdisreout", ir->nstdisreout);
1029 PR("orire-fc", ir->orires_fc);
1030 PR("orire-tau", ir->orires_tau);
1031 PR("nstorireout", ir->nstorireout);
1033 /* FREE ENERGY VARIABLES */
1034 PS("free-energy", enumValueToString(ir->efep));
1035 if (ir->efep != FreeEnergyPerturbationType::No || ir->bSimTemp)
1037 pr_fepvals(fp, indent, ir->fepvals.get(), bMDPformat);
1041 pr_expandedvals(fp, indent, ir->expandedvals.get(), ir->fepvals->n_lambda);
1044 /* NON-equilibrium MD stuff */
1045 PR("cos-acceleration", ir->cos_accel);
1046 pr_matrix(fp, indent, "deform", ir->deform, bMDPformat);
1048 /* SIMULATED TEMPERING */
1049 PS("simulated-tempering", EBOOL(ir->bSimTemp));
1052 pr_simtempvals(fp, indent, ir->simtempvals.get(), ir->fepvals->n_lambda);
1055 /* ION/WATER SWAPPING FOR COMPUTATIONAL ELECTROPHYSIOLOGY */
1056 PS("swapcoords", enumValueToString(ir->eSwapCoords));
1057 if (ir->eSwapCoords != SwapType::No)
1059 pr_swap(fp, indent, ir->swap);
1062 /* USER-DEFINED THINGIES */
1063 PI("userint1", ir->userint1);
1064 PI("userint2", ir->userint2);
1065 PI("userint3", ir->userint3);
1066 PI("userint4", ir->userint4);
1067 PR("userreal1", ir->userreal1);
1068 PR("userreal2", ir->userreal2);
1069 PR("userreal3", ir->userreal3);
1070 PR("userreal4", ir->userreal4);
1074 gmx::TextWriter writer(fp);
1075 writer.wrapperSettings().setIndent(indent);
1076 gmx::dumpKeyValueTree(&writer, *ir->params);
1079 pr_grp_opts(fp, indent, "grpopts", &(ir->opts), bMDPformat);
1086 static void cmp_grpopts(FILE* fp, const t_grpopts* opt1, const t_grpopts* opt2, real ftol, real abstol)
1089 char buf1[256], buf2[256];
1091 cmp_int(fp, "inputrec->grpopts.ngtc", -1, opt1->ngtc, opt2->ngtc);
1092 cmp_int(fp, "inputrec->grpopts.ngacc", -1, opt1->ngacc, opt2->ngacc);
1093 cmp_int(fp, "inputrec->grpopts.ngfrz", -1, opt1->ngfrz, opt2->ngfrz);
1094 cmp_int(fp, "inputrec->grpopts.ngener", -1, opt1->ngener, opt2->ngener);
1095 for (i = 0; (i < std::min(opt1->ngtc, opt2->ngtc)); i++)
1097 cmp_real(fp, "inputrec->grpopts.nrdf", i, opt1->nrdf[i], opt2->nrdf[i], ftol, abstol);
1098 cmp_real(fp, "inputrec->grpopts.ref_t", i, opt1->ref_t[i], opt2->ref_t[i], ftol, abstol);
1099 cmp_real(fp, "inputrec->grpopts.tau_t", i, opt1->tau_t[i], opt2->tau_t[i], ftol, abstol);
1100 cmpEnum(fp, "inputrec->grpopts.annealing", opt1->annealing[i], opt2->annealing[i]);
1101 cmp_int(fp, "inputrec->grpopts.anneal_npoints", i, opt1->anneal_npoints[i], opt2->anneal_npoints[i]);
1102 if (opt1->anneal_npoints[i] == opt2->anneal_npoints[i])
1104 sprintf(buf1, "inputrec->grpopts.anneal_time[%d]", i);
1105 sprintf(buf2, "inputrec->grpopts.anneal_temp[%d]", i);
1106 for (j = 0; j < opt1->anneal_npoints[i]; j++)
1108 cmp_real(fp, buf1, j, opt1->anneal_time[i][j], opt2->anneal_time[i][j], ftol, abstol);
1109 cmp_real(fp, buf2, j, opt1->anneal_temp[i][j], opt2->anneal_temp[i][j], ftol, abstol);
1113 if (opt1->ngener == opt2->ngener)
1115 for (i = 0; i < opt1->ngener; i++)
1117 for (j = i; j < opt1->ngener; j++)
1119 sprintf(buf1, "inputrec->grpopts.egp_flags[%d]", i);
1120 cmp_int(fp, buf1, j, opt1->egp_flags[opt1->ngener * i + j], opt2->egp_flags[opt1->ngener * i + j]);
1124 for (i = 0; (i < std::min(opt1->ngacc, opt2->ngacc)); i++)
1126 cmp_rvec(fp, "inputrec->grpopts.acceleration", i, opt1->acceleration[i], opt2->acceleration[i], ftol, abstol);
1128 for (i = 0; (i < std::min(opt1->ngfrz, opt2->ngfrz)); i++)
1130 cmp_ivec(fp, "inputrec->grpopts.nFreeze", i, opt1->nFreeze[i], opt2->nFreeze[i]);
1134 static void cmp_pull(FILE* fp)
1137 "WARNING: Both files use COM pulling, but comparing of the pull struct is not "
1138 "implemented (yet). The pull parameters could be the same or different.\n");
1141 static void cmp_awhDimParams(FILE* fp,
1142 const gmx::AwhDimParams& dimp1,
1143 const gmx::AwhDimParams& dimp2,
1148 /* Note that we have double index here, but the compare functions only
1149 * support one index, so here we only print the dim index and not the bias.
1152 "inputrec.awhParams->bias?->dim->coord_index",
1154 dimp1.coordinateIndex(),
1155 dimp2.coordinateIndex());
1156 cmp_double(fp, "inputrec->awhParams->bias?->dim->period", dimIndex, dimp1.period(), dimp2.period(), ftol, abstol);
1158 "inputrec->awhParams->bias?->dim->diffusion",
1164 cmp_double(fp, "inputrec->awhParams->bias?->dim->origin", dimIndex, dimp1.origin(), dimp2.origin(), ftol, abstol);
1165 cmp_double(fp, "inputrec->awhParams->bias?->dim->end", dimIndex, dimp1.end(), dimp2.end(), ftol, abstol);
1167 "inputrec->awhParams->bias?->dim->coord_value_init",
1169 dimp1.initialCoordinate(),
1170 dimp2.initialCoordinate(),
1174 "inputrec->awhParams->bias?->dim->coverDiameter",
1176 dimp1.coverDiameter(),
1177 dimp2.coverDiameter(),
1182 static void cmp_awhBiasParams(FILE* fp,
1183 const gmx::AwhBiasParams& bias1,
1184 const gmx::AwhBiasParams& bias2,
1189 cmp_int(fp, "inputrec->awhParams->ndim", biasIndex, bias1.ndim(), bias2.ndim());
1190 cmpEnum<gmx::AwhTargetType>(
1191 fp, "inputrec->awhParams->biaseTarget", bias1.targetDistribution(), bias2.targetDistribution());
1193 "inputrec->awhParams->biastargetBetaScaling",
1195 bias1.targetBetaScaling(),
1196 bias2.targetBetaScaling(),
1200 "inputrec->awhParams->biastargetCutoff",
1202 bias1.targetCutoff(),
1203 bias2.targetCutoff(),
1206 cmpEnum<gmx::AwhHistogramGrowthType>(
1207 fp, "inputrec->awhParams->biaseGrowth", bias1.growthType(), bias2.growthType());
1208 cmp_bool(fp, "inputrec->awhParams->biasbUserData", biasIndex, bias1.userPMFEstimate(), bias2.userPMFEstimate());
1210 "inputrec->awhParams->biaserror_initial",
1212 bias1.initialErrorEstimate(),
1213 bias2.initialErrorEstimate(),
1216 cmp_int(fp, "inputrec->awhParams->biasShareGroup", biasIndex, bias1.shareGroup(), bias2.shareGroup());
1218 const auto dimParams1 = bias1.dimParams();
1219 const auto dimParams2 = bias2.dimParams();
1220 for (int dim = 0; dim < std::min(bias1.ndim(), bias2.ndim()); dim++)
1222 cmp_awhDimParams(fp, dimParams1[dim], dimParams2[dim], dim, ftol, abstol);
1226 static void cmp_awhParams(FILE* fp, const gmx::AwhParams& awh1, const gmx::AwhParams& awh2, real ftol, real abstol)
1228 cmp_int(fp, "inputrec->awhParams->nbias", -1, awh1.numBias(), awh2.numBias());
1229 cmp_int64(fp, "inputrec->awhParams->seed", awh1.seed(), awh2.seed());
1230 cmp_int(fp, "inputrec->awhParams->nstout", -1, awh1.nstout(), awh2.nstout());
1231 cmp_int(fp, "inputrec->awhParams->nstsample_coord", -1, awh1.nstSampleCoord(), awh2.nstSampleCoord());
1233 "inputrec->awhParams->nsamples_update_free_energy",
1235 awh1.numSamplesUpdateFreeEnergy(),
1236 awh2.numSamplesUpdateFreeEnergy());
1237 cmpEnum<gmx::AwhPotentialType>(
1238 fp, "inputrec->awhParams->ePotential", awh1.potential(), awh2.potential());
1239 cmp_bool(fp, "inputrec->awhParams->shareBiasMultisim", -1, awh1.shareBiasMultisim(), awh2.shareBiasMultisim());
1241 if (awh1.numBias() == awh2.numBias())
1243 const auto awhBiasParams1 = awh1.awhBiasParams();
1244 const auto awhBiasParams2 = awh2.awhBiasParams();
1245 for (int bias = 0; bias < awh1.numBias(); bias++)
1247 cmp_awhBiasParams(fp, awhBiasParams1[bias], awhBiasParams2[bias], bias, ftol, abstol);
1252 static void cmp_simtempvals(FILE* fp,
1253 const t_simtemp* simtemp1,
1254 const t_simtemp* simtemp2,
1260 cmpEnum(fp, "inputrec->simtempvals->eSimTempScale", simtemp1->eSimTempScale, simtemp2->eSimTempScale);
1261 cmp_real(fp, "inputrec->simtempvals->simtemp_high", -1, simtemp1->simtemp_high, simtemp2->simtemp_high, ftol, abstol);
1262 cmp_real(fp, "inputrec->simtempvals->simtemp_low", -1, simtemp1->simtemp_low, simtemp2->simtemp_low, ftol, abstol);
1263 for (i = 0; i < n_lambda; i++)
1266 "inputrec->simtempvals->temperatures",
1268 simtemp1->temperatures[i],
1269 simtemp2->temperatures[i],
1275 static void cmp_expandedvals(FILE* fp,
1276 const t_expanded* expand1,
1277 const t_expanded* expand2,
1284 cmp_bool(fp, "inputrec->fepvals->bInit_weights", -1, expand1->bInit_weights, expand2->bInit_weights);
1285 cmp_bool(fp, "inputrec->fepvals->bWLoneovert", -1, expand1->bWLoneovert, expand2->bWLoneovert);
1287 for (i = 0; i < n_lambda; i++)
1290 "inputrec->expandedvals->init_lambda_weights",
1292 expand1->init_lambda_weights[i],
1293 expand2->init_lambda_weights[i],
1298 cmpEnum(fp, "inputrec->expandedvals->lambda-stats", expand1->elamstats, expand2->elamstats);
1299 cmpEnum(fp, "inputrec->expandedvals->lambda-mc-move", expand1->elmcmove, expand2->elmcmove);
1300 cmp_int(fp, "inputrec->expandedvals->lmc-repeats", -1, expand1->lmc_repeats, expand2->lmc_repeats);
1301 cmp_int(fp, "inputrec->expandedvals->lmc-gibbsdelta", -1, expand1->gibbsdeltalam, expand2->gibbsdeltalam);
1302 cmp_int(fp, "inputrec->expandedvals->lmc-forced-nstart", -1, expand1->lmc_forced_nstart, expand2->lmc_forced_nstart);
1303 cmpEnum(fp, "inputrec->expandedvals->lambda-weights-equil", expand1->elmceq, expand2->elmceq);
1305 "inputrec->expandedvals->,weight-equil-number-all-lambda",
1307 expand1->equil_n_at_lam,
1308 expand2->equil_n_at_lam);
1309 cmp_int(fp, "inputrec->expandedvals->weight-equil-number-samples", -1, expand1->equil_samples, expand2->equil_samples);
1310 cmp_int(fp, "inputrec->expandedvals->weight-equil-number-steps", -1, expand1->equil_steps, expand2->equil_steps);
1312 "inputrec->expandedvals->weight-equil-wl-delta",
1314 expand1->equil_wl_delta,
1315 expand2->equil_wl_delta,
1319 "inputrec->expandedvals->weight-equil-count-ratio",
1321 expand1->equil_ratio,
1322 expand2->equil_ratio,
1326 "inputrec->expandedvals->symmetrized-transition-matrix",
1328 expand1->bSymmetrizedTMatrix,
1329 expand2->bSymmetrizedTMatrix);
1330 cmp_int(fp, "inputrec->expandedvals->nstTij", -1, expand1->nstTij, expand2->nstTij);
1331 cmp_int(fp, "inputrec->expandedvals->mininum-var-min", -1, expand1->minvarmin, expand2->minvarmin); /*default is reasonable */
1332 cmp_int(fp, "inputrec->expandedvals->weight-c-range", -1, expand1->c_range, expand2->c_range); /* default is just C=0 */
1333 cmp_real(fp, "inputrec->expandedvals->wl-scale", -1, expand1->wl_scale, expand2->wl_scale, ftol, abstol);
1334 cmp_real(fp, "inputrec->expandedvals->init-wl-delta", -1, expand1->init_wl_delta, expand2->init_wl_delta, ftol, abstol);
1335 cmp_real(fp, "inputrec->expandedvals->wl-ratio", -1, expand1->wl_ratio, expand2->wl_ratio, ftol, abstol);
1336 cmp_int(fp, "inputrec->expandedvals->nstexpanded", -1, expand1->nstexpanded, expand2->nstexpanded);
1337 cmp_int(fp, "inputrec->expandedvals->lmc-seed", -1, expand1->lmc_seed, expand2->lmc_seed);
1338 cmp_real(fp, "inputrec->expandedvals->mc-temperature", -1, expand1->mc_temp, expand2->mc_temp, ftol, abstol);
1341 static void cmp_fepvals(FILE* fp, const t_lambda* fep1, const t_lambda* fep2, real ftol, real abstol)
1344 cmp_int(fp, "inputrec->nstdhdl", -1, fep1->nstdhdl, fep2->nstdhdl);
1345 cmp_double(fp, "inputrec->fepvals->init_fep_state", -1, fep1->init_fep_state, fep2->init_fep_state, ftol, abstol);
1346 cmp_double(fp, "inputrec->fepvals->delta_lambda", -1, fep1->delta_lambda, fep2->delta_lambda, ftol, abstol);
1347 cmp_int(fp, "inputrec->fepvals->n_lambda", -1, fep1->n_lambda, fep2->n_lambda);
1348 for (i = 0; i < static_cast<int>(FreeEnergyPerturbationCouplingType::Count); i++)
1350 for (j = 0; j < std::min(fep1->n_lambda, fep2->n_lambda); j++)
1353 "inputrec->fepvals->all_lambda",
1355 fep1->all_lambda[i][j],
1356 fep2->all_lambda[i][j],
1361 cmp_int(fp, "inputrec->fepvals->lambda_neighbors", 1, fep1->lambda_neighbors, fep2->lambda_neighbors);
1362 cmp_real(fp, "inputrec->fepvals->sc_alpha", -1, fep1->sc_alpha, fep2->sc_alpha, ftol, abstol);
1363 cmp_int(fp, "inputrec->fepvals->sc_power", -1, fep1->sc_power, fep2->sc_power);
1364 cmp_real(fp, "inputrec->fepvals->sc_r_power", -1, fep1->sc_r_power, fep2->sc_r_power, ftol, abstol);
1365 cmp_real(fp, "inputrec->fepvals->sc_sigma", -1, fep1->sc_sigma, fep2->sc_sigma, ftol, abstol);
1366 cmpEnum(fp, "inputrec->fepvals->edHdLPrintEnergy", fep1->edHdLPrintEnergy, fep1->edHdLPrintEnergy);
1367 cmp_bool(fp, "inputrec->fepvals->bScCoul", -1, fep1->bScCoul, fep1->bScCoul);
1368 cmpEnum(fp, "inputrec->separate_dhdl_file", fep1->separate_dhdl_file, fep2->separate_dhdl_file);
1369 cmpEnum(fp, "inputrec->dhdl_derivatives", fep1->dhdl_derivatives, fep2->dhdl_derivatives);
1370 cmp_int(fp, "inputrec->dh_hist_size", -1, fep1->dh_hist_size, fep2->dh_hist_size);
1371 cmp_double(fp, "inputrec->dh_hist_spacing", -1, fep1->dh_hist_spacing, fep2->dh_hist_spacing, ftol, abstol);
1372 cmpEnum(fp, "inputrec->fepvals->softcoreFunction", fep1->softcoreFunction, fep2->softcoreFunction);
1374 "inputrec->fepvals->scGapsysScaleLinpointLJ",
1376 fep1->scGapsysScaleLinpointLJ,
1377 fep2->scGapsysScaleLinpointLJ,
1381 "inputrec->fepvals->scGapsysScaleLinpointQ",
1383 fep1->scGapsysScaleLinpointQ,
1384 fep2->scGapsysScaleLinpointQ,
1387 cmp_real(fp, "inputrec->fepvals->scGapsysSigmaLJ", -1, fep1->scGapsysSigmaLJ, fep2->scGapsysSigmaLJ, ftol, abstol);
1390 void cmp_inputrec(FILE* fp, const t_inputrec* ir1, const t_inputrec* ir2, real ftol, real abstol)
1392 fprintf(fp, "comparing inputrec\n");
1394 /* gcc 2.96 doesnt like these defines at all, but issues a huge list
1395 * of warnings. Maybe it will change in future versions, but for the
1396 * moment I've spelled them out instead. /EL 000820
1397 * #define CIB(s) cmp_int(fp,"inputrec->"#s,0,ir1->##s,ir2->##s)
1398 * #define CII(s) cmp_int(fp,"inputrec->"#s,0,ir1->##s,ir2->##s)
1399 * #define CIR(s) cmp_real(fp,"inputrec->"#s,0,ir1->##s,ir2->##s,ftol)
1401 cmpEnum(fp, "inputrec->eI", ir1->eI, ir2->eI);
1402 cmp_int64(fp, "inputrec->nsteps", ir1->nsteps, ir2->nsteps);
1403 cmp_int64(fp, "inputrec->init_step", ir1->init_step, ir2->init_step);
1404 cmp_int(fp, "inputrec->simulation_part", -1, ir1->simulation_part, ir2->simulation_part);
1405 cmp_int(fp, "inputrec->mts", -1, static_cast<int>(ir1->useMts), static_cast<int>(ir2->useMts));
1406 if (ir1->useMts && ir2->useMts)
1408 cmp_int(fp, "inputrec->mts-levels", -1, ir1->mtsLevels.size(), ir2->mtsLevels.size());
1410 "inputrec->mts-level2-forces",
1412 ir1->mtsLevels[1].forceGroups.to_ulong(),
1413 ir2->mtsLevels[1].forceGroups.to_ulong());
1415 "inputrec->mts-level2-factor",
1417 ir1->mtsLevels[1].stepFactor,
1418 ir2->mtsLevels[1].stepFactor);
1420 cmp_int(fp, "inputrec->pbcType", -1, static_cast<int>(ir1->pbcType), static_cast<int>(ir2->pbcType));
1421 cmp_bool(fp, "inputrec->bPeriodicMols", -1, ir1->bPeriodicMols, ir2->bPeriodicMols);
1422 cmpEnum(fp, "inputrec->cutoff_scheme", ir1->cutoff_scheme, ir2->cutoff_scheme);
1423 cmp_int(fp, "inputrec->nstlist", -1, ir1->nstlist, ir2->nstlist);
1424 cmp_int(fp, "inputrec->nstcomm", -1, ir1->nstcomm, ir2->nstcomm);
1425 cmpEnum(fp, "inputrec->comm_mode", ir1->comm_mode, ir2->comm_mode);
1426 cmp_int(fp, "inputrec->nstlog", -1, ir1->nstlog, ir2->nstlog);
1427 cmp_int(fp, "inputrec->nstxout", -1, ir1->nstxout, ir2->nstxout);
1428 cmp_int(fp, "inputrec->nstvout", -1, ir1->nstvout, ir2->nstvout);
1429 cmp_int(fp, "inputrec->nstfout", -1, ir1->nstfout, ir2->nstfout);
1430 cmp_int(fp, "inputrec->nstcalcenergy", -1, ir1->nstcalcenergy, ir2->nstcalcenergy);
1431 cmp_int(fp, "inputrec->nstenergy", -1, ir1->nstenergy, ir2->nstenergy);
1432 cmp_int(fp, "inputrec->nstxout_compressed", -1, ir1->nstxout_compressed, ir2->nstxout_compressed);
1433 cmp_double(fp, "inputrec->init_t", -1, ir1->init_t, ir2->init_t, ftol, abstol);
1434 cmp_double(fp, "inputrec->delta_t", -1, ir1->delta_t, ir2->delta_t, ftol, abstol);
1436 "inputrec->x_compression_precision",
1438 ir1->x_compression_precision,
1439 ir2->x_compression_precision,
1442 cmp_real(fp, "inputrec->fourierspacing", -1, ir1->fourier_spacing, ir2->fourier_spacing, ftol, abstol);
1443 cmp_int(fp, "inputrec->nkx", -1, ir1->nkx, ir2->nkx);
1444 cmp_int(fp, "inputrec->nky", -1, ir1->nky, ir2->nky);
1445 cmp_int(fp, "inputrec->nkz", -1, ir1->nkz, ir2->nkz);
1446 cmp_int(fp, "inputrec->pme_order", -1, ir1->pme_order, ir2->pme_order);
1447 cmp_real(fp, "inputrec->ewald_rtol", -1, ir1->ewald_rtol, ir2->ewald_rtol, ftol, abstol);
1448 cmpEnum(fp, "inputrec->ewald_geometry", ir1->ewald_geometry, ir2->ewald_geometry);
1449 cmp_real(fp, "inputrec->epsilon_surface", -1, ir1->epsilon_surface, ir2->epsilon_surface, ftol, abstol);
1451 "inputrec->bContinuation",
1453 static_cast<int>(ir1->bContinuation),
1454 static_cast<int>(ir2->bContinuation));
1455 cmp_int(fp, "inputrec->bShakeSOR", -1, static_cast<int>(ir1->bShakeSOR), static_cast<int>(ir2->bShakeSOR));
1456 cmpEnum(fp, "inputrec->etc", ir1->etc, ir2->etc);
1458 "inputrec->bPrintNHChains",
1460 static_cast<int>(ir1->bPrintNHChains),
1461 static_cast<int>(ir2->bPrintNHChains));
1462 cmpEnum(fp, "inputrec->epc", ir1->epc, ir2->epc);
1463 cmpEnum(fp, "inputrec->epct", ir1->epct, ir2->epct);
1464 cmp_real(fp, "inputrec->tau_p", -1, ir1->tau_p, ir2->tau_p, ftol, abstol);
1465 cmp_rvec(fp, "inputrec->ref_p(x)", -1, ir1->ref_p[XX], ir2->ref_p[XX], ftol, abstol);
1466 cmp_rvec(fp, "inputrec->ref_p(y)", -1, ir1->ref_p[YY], ir2->ref_p[YY], ftol, abstol);
1467 cmp_rvec(fp, "inputrec->ref_p(z)", -1, ir1->ref_p[ZZ], ir2->ref_p[ZZ], ftol, abstol);
1468 cmp_rvec(fp, "inputrec->compress(x)", -1, ir1->compress[XX], ir2->compress[XX], ftol, abstol);
1469 cmp_rvec(fp, "inputrec->compress(y)", -1, ir1->compress[YY], ir2->compress[YY], ftol, abstol);
1470 cmp_rvec(fp, "inputrec->compress(z)", -1, ir1->compress[ZZ], ir2->compress[ZZ], ftol, abstol);
1471 cmpEnum(fp, "refcoord_scaling", ir1->refcoord_scaling, ir2->refcoord_scaling);
1472 cmp_rvec(fp, "inputrec->posres_com", -1, ir1->posres_com, ir2->posres_com, ftol, abstol);
1473 cmp_rvec(fp, "inputrec->posres_comB", -1, ir1->posres_comB, ir2->posres_comB, ftol, abstol);
1474 cmp_real(fp, "inputrec->verletbuf_tol", -1, ir1->verletbuf_tol, ir2->verletbuf_tol, ftol, abstol);
1475 cmp_real(fp, "inputrec->rlist", -1, ir1->rlist, ir2->rlist, ftol, abstol);
1476 cmp_real(fp, "inputrec->rtpi", -1, ir1->rtpi, ir2->rtpi, ftol, abstol);
1477 cmpEnum(fp, "inputrec->coulombtype", ir1->coulombtype, ir2->coulombtype);
1478 cmpEnum(fp, "inputrec->coulomb_modifier", ir1->coulomb_modifier, ir2->coulomb_modifier);
1479 cmp_real(fp, "inputrec->rcoulomb_switch", -1, ir1->rcoulomb_switch, ir2->rcoulomb_switch, ftol, abstol);
1480 cmp_real(fp, "inputrec->rcoulomb", -1, ir1->rcoulomb, ir2->rcoulomb, ftol, abstol);
1481 cmpEnum(fp, "inputrec->vdwtype", ir1->vdwtype, ir2->vdwtype);
1482 cmpEnum(fp, "inputrec->vdw_modifier", ir1->vdw_modifier, ir2->vdw_modifier);
1483 cmp_real(fp, "inputrec->rvdw_switch", -1, ir1->rvdw_switch, ir2->rvdw_switch, ftol, abstol);
1484 cmp_real(fp, "inputrec->rvdw", -1, ir1->rvdw, ir2->rvdw, ftol, abstol);
1485 cmp_real(fp, "inputrec->epsilon_r", -1, ir1->epsilon_r, ir2->epsilon_r, ftol, abstol);
1486 cmp_real(fp, "inputrec->epsilon_rf", -1, ir1->epsilon_rf, ir2->epsilon_rf, ftol, abstol);
1487 cmp_real(fp, "inputrec->tabext", -1, ir1->tabext, ir2->tabext, ftol, abstol);
1489 cmpEnum(fp, "inputrec->eDispCorr", ir1->eDispCorr, ir2->eDispCorr);
1490 cmp_real(fp, "inputrec->shake_tol", -1, ir1->shake_tol, ir2->shake_tol, ftol, abstol);
1491 cmpEnum(fp, "inputrec->efep", ir1->efep, ir2->efep);
1492 cmp_fepvals(fp, ir1->fepvals.get(), ir2->fepvals.get(), ftol, abstol);
1493 cmp_int(fp, "inputrec->bSimTemp", -1, static_cast<int>(ir1->bSimTemp), static_cast<int>(ir2->bSimTemp));
1494 if ((ir1->bSimTemp == ir2->bSimTemp) && (ir1->bSimTemp))
1497 ir1->simtempvals.get(),
1498 ir2->simtempvals.get(),
1499 std::min(ir1->fepvals->n_lambda, ir2->fepvals->n_lambda),
1503 cmp_int(fp, "inputrec->bExpanded", -1, static_cast<int>(ir1->bExpanded), static_cast<int>(ir2->bExpanded));
1504 if ((ir1->bExpanded == ir2->bExpanded) && (ir1->bExpanded))
1506 cmp_expandedvals(fp,
1507 ir1->expandedvals.get(),
1508 ir2->expandedvals.get(),
1509 std::min(ir1->fepvals->n_lambda, ir2->fepvals->n_lambda),
1513 cmp_int(fp, "inputrec->nwall", -1, ir1->nwall, ir2->nwall);
1514 cmpEnum(fp, "inputrec->wall_type", ir1->wall_type, ir2->wall_type);
1515 cmp_int(fp, "inputrec->wall_atomtype[0]", -1, ir1->wall_atomtype[0], ir2->wall_atomtype[0]);
1516 cmp_int(fp, "inputrec->wall_atomtype[1]", -1, ir1->wall_atomtype[1], ir2->wall_atomtype[1]);
1517 cmp_real(fp, "inputrec->wall_density[0]", -1, ir1->wall_density[0], ir2->wall_density[0], ftol, abstol);
1518 cmp_real(fp, "inputrec->wall_density[1]", -1, ir1->wall_density[1], ir2->wall_density[1], ftol, abstol);
1519 cmp_real(fp, "inputrec->wall_ewald_zfac", -1, ir1->wall_ewald_zfac, ir2->wall_ewald_zfac, ftol, abstol);
1521 cmp_bool(fp, "inputrec->bPull", -1, ir1->bPull, ir2->bPull);
1522 if (ir1->bPull && ir2->bPull)
1527 cmp_bool(fp, "inputrec->bDoAwh", -1, ir1->bDoAwh, ir2->bDoAwh);
1528 if (ir1->bDoAwh && ir2->bDoAwh)
1530 cmp_awhParams(fp, *ir1->awhParams, *ir2->awhParams, ftol, abstol);
1533 cmpEnum(fp, "inputrec->eDisre", ir1->eDisre, ir2->eDisre);
1534 cmp_real(fp, "inputrec->dr_fc", -1, ir1->dr_fc, ir2->dr_fc, ftol, abstol);
1535 cmpEnum(fp, "inputrec->eDisreWeighting", ir1->eDisreWeighting, ir2->eDisreWeighting);
1536 cmp_int(fp, "inputrec->bDisreMixed", -1, static_cast<int>(ir1->bDisreMixed), static_cast<int>(ir2->bDisreMixed));
1537 cmp_int(fp, "inputrec->nstdisreout", -1, ir1->nstdisreout, ir2->nstdisreout);
1538 cmp_real(fp, "inputrec->dr_tau", -1, ir1->dr_tau, ir2->dr_tau, ftol, abstol);
1539 cmp_real(fp, "inputrec->orires_fc", -1, ir1->orires_fc, ir2->orires_fc, ftol, abstol);
1540 cmp_real(fp, "inputrec->orires_tau", -1, ir1->orires_tau, ir2->orires_tau, ftol, abstol);
1541 cmp_int(fp, "inputrec->nstorireout", -1, ir1->nstorireout, ir2->nstorireout);
1542 cmp_real(fp, "inputrec->em_stepsize", -1, ir1->em_stepsize, ir2->em_stepsize, ftol, abstol);
1543 cmp_real(fp, "inputrec->em_tol", -1, ir1->em_tol, ir2->em_tol, ftol, abstol);
1544 cmp_int(fp, "inputrec->niter", -1, ir1->niter, ir2->niter);
1545 cmp_real(fp, "inputrec->fc_stepsize", -1, ir1->fc_stepsize, ir2->fc_stepsize, ftol, abstol);
1546 cmp_int(fp, "inputrec->nstcgsteep", -1, ir1->nstcgsteep, ir2->nstcgsteep);
1547 cmp_int(fp, "inputrec->nbfgscorr", 0, ir1->nbfgscorr, ir2->nbfgscorr);
1548 cmpEnum(fp, "inputrec->eConstrAlg", ir1->eConstrAlg, ir2->eConstrAlg);
1549 cmp_int(fp, "inputrec->nProjOrder", -1, ir1->nProjOrder, ir2->nProjOrder);
1550 cmp_real(fp, "inputrec->LincsWarnAngle", -1, ir1->LincsWarnAngle, ir2->LincsWarnAngle, ftol, abstol);
1551 cmp_int(fp, "inputrec->nLincsIter", -1, ir1->nLincsIter, ir2->nLincsIter);
1552 cmp_real(fp, "inputrec->bd_fric", -1, ir1->bd_fric, ir2->bd_fric, ftol, abstol);
1553 cmp_int64(fp, "inputrec->ld_seed", ir1->ld_seed, ir2->ld_seed);
1554 cmp_real(fp, "inputrec->cos_accel", -1, ir1->cos_accel, ir2->cos_accel, ftol, abstol);
1555 cmp_rvec(fp, "inputrec->deform(a)", -1, ir1->deform[XX], ir2->deform[XX], ftol, abstol);
1556 cmp_rvec(fp, "inputrec->deform(b)", -1, ir1->deform[YY], ir2->deform[YY], ftol, abstol);
1557 cmp_rvec(fp, "inputrec->deform(c)", -1, ir1->deform[ZZ], ir2->deform[ZZ], ftol, abstol);
1560 cmp_int(fp, "inputrec->userint1", -1, ir1->userint1, ir2->userint1);
1561 cmp_int(fp, "inputrec->userint2", -1, ir1->userint2, ir2->userint2);
1562 cmp_int(fp, "inputrec->userint3", -1, ir1->userint3, ir2->userint3);
1563 cmp_int(fp, "inputrec->userint4", -1, ir1->userint4, ir2->userint4);
1564 cmp_real(fp, "inputrec->userreal1", -1, ir1->userreal1, ir2->userreal1, ftol, abstol);
1565 cmp_real(fp, "inputrec->userreal2", -1, ir1->userreal2, ir2->userreal2, ftol, abstol);
1566 cmp_real(fp, "inputrec->userreal3", -1, ir1->userreal3, ir2->userreal3, ftol, abstol);
1567 cmp_real(fp, "inputrec->userreal4", -1, ir1->userreal4, ir2->userreal4, ftol, abstol);
1568 cmp_grpopts(fp, &(ir1->opts), &(ir2->opts), ftol, abstol);
1569 gmx::TextWriter writer(fp);
1570 gmx::compareKeyValueTrees(&writer, *ir1->params, *ir2->params, ftol, abstol);
1573 void comp_pull_AB(FILE* fp, const pull_params_t& pull, real ftol, real abstol)
1575 for (int i = 0; i < pull.ncoord; i++)
1577 fprintf(fp, "comparing pull coord %d\n", i);
1578 cmp_real(fp, "pull-coord->k", -1, pull.coord[i].k, pull.coord[i].kB, ftol, abstol);
1582 gmx_bool inputrecDeform(const t_inputrec* ir)
1584 return (ir->deform[XX][XX] != 0 || ir->deform[YY][YY] != 0 || ir->deform[ZZ][ZZ] != 0
1585 || ir->deform[YY][XX] != 0 || ir->deform[ZZ][XX] != 0 || ir->deform[ZZ][YY] != 0);
1588 gmx_bool inputrecDynamicBox(const t_inputrec* ir)
1590 return (ir->epc != PressureCoupling::No || ir->eI == IntegrationAlgorithm::TPI || inputrecDeform(ir));
1593 gmx_bool inputrecPreserveShape(const t_inputrec* ir)
1595 return (ir->epc != PressureCoupling::No && ir->deform[XX][XX] == 0
1596 && (ir->epct == PressureCouplingType::Isotropic
1597 || ir->epct == PressureCouplingType::SemiIsotropic));
1600 gmx_bool inputrecNeedMutot(const t_inputrec* ir)
1602 return ((ir->coulombtype == CoulombInteractionType::Ewald || EEL_PME(ir->coulombtype))
1603 && (ir->ewald_geometry == EwaldGeometry::ThreeDC || ir->epsilon_surface != 0));
1606 gmx_bool inputrecExclForces(const t_inputrec* ir)
1608 return (EEL_FULL(ir->coulombtype) || (EEL_RF(ir->coulombtype)));
1611 gmx_bool inputrecNptTrotter(const t_inputrec* ir)
1613 return (((ir->eI == IntegrationAlgorithm::VV) || (ir->eI == IntegrationAlgorithm::VVAK))
1614 && (ir->epc == PressureCoupling::Mttk) && (ir->etc == TemperatureCoupling::NoseHoover));
1617 gmx_bool inputrecNvtTrotter(const t_inputrec* ir)
1619 return (((ir->eI == IntegrationAlgorithm::VV) || (ir->eI == IntegrationAlgorithm::VVAK))
1620 && (ir->epc != PressureCoupling::Mttk) && (ir->etc == TemperatureCoupling::NoseHoover));
1623 gmx_bool inputrecNphTrotter(const t_inputrec* ir)
1625 return (((ir->eI == IntegrationAlgorithm::VV) || (ir->eI == IntegrationAlgorithm::VVAK))
1626 && (ir->epc == PressureCoupling::Mttk) && (ir->etc != TemperatureCoupling::NoseHoover));
1629 bool inputrecPbcXY2Walls(const t_inputrec* ir)
1631 return (ir->pbcType == PbcType::XY && ir->nwall == 2);
1634 bool inputrecFrozenAtoms(const t_inputrec* ir)
1636 return ((ir->opts.nFreeze != nullptr)
1637 && (ir->opts.ngfrz > 1 || ir->opts.nFreeze[0][XX] != 0 || ir->opts.nFreeze[0][YY] != 0
1638 || ir->opts.nFreeze[0][ZZ] != 0));
1641 bool integratorHasConservedEnergyQuantity(const t_inputrec* ir)
1644 { // NOLINT bugprone-branch-clone
1645 // Energy minimization or stochastic integrator: no conservation
1648 else if (ir->etc == TemperatureCoupling::No && ir->epc == PressureCoupling::No)
1650 // The total energy is conserved, no additional conserved quanitity
1655 // Shear stress with Parrinello-Rahman is not supported (tedious)
1657 ((ir->epc == PressureCoupling::ParrinelloRahman || ir->epc == PressureCoupling::Mttk)
1658 && (ir->ref_p[YY][XX] != 0 || ir->ref_p[ZZ][XX] != 0 || ir->ref_p[ZZ][YY] != 0));
1660 return !ETC_ANDERSEN(ir->etc) && !shearWithPR;
1664 bool integratorHasReferenceTemperature(const t_inputrec* ir)
1666 return ((ir->etc != TemperatureCoupling::No) || EI_SD(ir->eI)
1667 || (ir->eI == IntegrationAlgorithm::BD) || EI_TPI(ir->eI));
1670 int inputrec2nboundeddim(const t_inputrec* ir)
1672 if (inputrecPbcXY2Walls(ir))
1678 return numPbcDimensions(ir->pbcType);
1682 int ndof_com(const t_inputrec* ir)
1686 switch (ir->pbcType)
1689 case PbcType::No: n = 3; break;
1690 case PbcType::XY: n = (ir->nwall == 0 ? 3 : 2); break;
1691 case PbcType::Screw: n = 1; break;
1692 default: gmx_incons("Unknown pbc in calc_nrdf");
1698 real maxReferenceTemperature(const t_inputrec& ir)
1700 if (EI_ENERGY_MINIMIZATION(ir.eI) || ir.eI == IntegrationAlgorithm::NM)
1705 if (EI_MD(ir.eI) && ir.etc == TemperatureCoupling::No)
1710 /* SD and BD also use ref_t and tau_t for setting the reference temperature.
1711 * TPI can be treated as MD, since it needs an ensemble temperature.
1714 real maxTemperature = 0;
1715 for (int i = 0; i < ir.opts.ngtc; i++)
1717 if (ir.opts.tau_t[i] >= 0)
1719 maxTemperature = std::max(maxTemperature, ir.opts.ref_t[i]);
1723 return maxTemperature;
1726 bool haveEwaldSurfaceContribution(const t_inputrec& ir)
1728 return EEL_PME_EWALD(ir.coulombtype)
1729 && (ir.ewald_geometry == EwaldGeometry::ThreeDC || ir.epsilon_surface != 0);
1732 bool haveFreeEnergyType(const t_inputrec& ir, const int fepType)
1734 for (int i = 0; i < ir.fepvals->n_lambda; i++)
1736 if (ir.fepvals->all_lambda[fepType][i] > 0)