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-2004, The GROMACS development team.
6 * Copyright (c) 2013,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.
38 #ifndef GMX_MDTYPES_INPUTREC_H
39 #define GMX_MDTYPES_INPUTREC_H
46 #include "gromacs/math/vectypes.h"
47 #include "gromacs/mdtypes/md_enums.h"
48 #include "gromacs/utility/basedefinitions.h"
49 #include "gromacs/utility/enumerationhelpers.h"
50 #include "gromacs/utility/real.h"
52 #define EGP_EXCL (1 << 0)
53 #define EGP_TABLE (1 << 1)
63 class KeyValueTreeObject;
71 //! Number of T-Coupl groups
73 //! Number of of Nose-Hoover chains per group
75 //! Number of Freeze groups
77 //! Number of Energy groups
79 //! Number of degrees of freedom in a group
81 //! Coupling temperature per group
83 //! No/simple/periodic simulated annealing for each group
84 SimulatedAnnealing* annealing;
85 //! Number of annealing time points per group
87 //! For each group: Time points
89 //! For each group: Temperature at these times. Final temp after all intervals is ref_t
93 //! Whether the group will be frozen in each direction
95 //! Exclusions/tables of energy group pairs
99 //! Number of QM groups
105 //! Simulated temperature scaling; linear or exponential
106 SimulatedTempering eSimTempScale;
107 //! The low temperature for simulated tempering
109 //! The high temperature for simulated tempering
111 //! The range of temperatures used for simulated tempering
117 //! The frequency for calculating dhdl
119 //! Fractional value of lambda (usually will use init_fep_state, this will only be for slow growth, and for legacy free energy code. Only has a valid value if positive)
121 //! The initial number of the state
123 //! Change of lambda per time step (fraction of (0.1)
125 //! Print no, total or potential energies in dhdl
126 FreeEnergyPrintEnergy edHdLPrintEnergy;
127 //! The number of foreign lambda points
129 //! The array of all lambda values
130 gmx::EnumerationArray<FreeEnergyPerturbationCouplingType, std::vector<double>> all_lambda;
131 //! The number of neighboring lambda states to calculate the energy for in up and down directions (-1 for all)
132 int lambda_neighbors;
133 //! The first lambda to calculate energies for
135 //! The last lambda +1 to calculate energies for
137 //! Free energy soft-core parameter
139 //! Lambda power for soft-core interactions
141 //! R power for soft-core interactions
143 //! Free energy soft-core sigma when c6 or c12=0
145 //! Free energy soft-core sigma for ?????
147 //! Use softcore for the coulomb portion as well (default FALSE)
149 //! Whether to print the dvdl term associated with this term; if it is not specified as separate, it is lumped with the FEP term
150 gmx::EnumerationArray<FreeEnergyPerturbationCouplingType, bool> separate_dvdl;
151 //! Whether to write a separate dhdl.xvg file note: NOT a gmx_bool, but an enum
152 SeparateDhdlFile separate_dhdl_file;
153 //! Whether to calculate+write dhdl derivatives note: NOT a gmx_bool, but an enum
154 DhDlDerivativeCalculation dhdl_derivatives;
155 //! The maximum table size for the dH histogram
157 //! The spacing for the dH histogram
158 double dh_hist_spacing;
163 //! The frequency of expanded ensemble state changes
165 //! Which type of move updating do we use for lambda monte carlo (or no for none)
166 LambdaWeightCalculation elamstats;
167 //! What move set will be we using for state space moves
168 LambdaMoveCalculation elmcmove;
169 //! The method we use to decide of we have equilibrated the weights
170 LambdaWeightWillReachEquilibrium elmceq;
171 //! The minumum number of samples at each lambda for deciding whether we have reached a minimum
173 //! Wang-Landau delta at which we stop equilibrating weights
175 //! Use the ratio of weights (ratio of minimum to maximum) to decide when to stop equilibrating
177 //! After equil_steps steps we stop equilibrating the weights
179 //! After equil_samples total samples (steps/nstfep), we stop equilibrating the weights
181 //! Random number seed for lambda mc switches
183 //! Whether to use minumum variance weighting
185 //! The number of samples needed before kicking into minvar routine
187 //! The offset for the variance in MinVar
189 //! Range of cvalues used for BAR
191 //! Whether to print symmetrized matrices
192 gmx_bool bSymmetrizedTMatrix;
193 //! How frequently to print the transition matrices
195 //! Number of repetitions in the MC lambda jumps MRS -- VERIFY THIS
197 //! Minimum number of samples for each state before free sampling MRS -- VERIFY THIS!
198 int lmc_forced_nstart;
199 //! Distance in lambda space for the gibbs interval
201 //! Scaling factor for Wang-Landau
203 //! Ratio between largest and smallest number for freezing the weights
205 //! Starting delta for Wang-Landau
207 //! Use one over t convergence for Wang-Landau when the delta get sufficiently small
208 gmx_bool bWLoneovert;
209 //! Did we initialize the weights? TODO: REMOVE FOR 5.0, no longer needed with new logic
210 gmx_bool bInit_weights;
211 //! To override the main temperature, or define it if it's not defined
213 //! User-specified initial weights to start with
214 std::vector<real> init_lambda_weights;
219 //! Rotation type for this group
220 EnforcedRotationGroupType eType;
221 //! Use mass-weighed positions?
223 //! Number of atoms in the group
225 //! The global atoms numbers
227 //! The reference positions
229 //! The normalized rotation vector
231 //! Rate of rotation (degree/ps)
233 //! Force constant (kJ/(mol nm^2)
235 //! Pivot point of rotation axis (nm)
237 //! Type of fit to determine actual group angle
238 RotationGroupFitting eFittype;
239 //! Number of angles around the reference angle for which the rotation potential is also evaluated (for fit type 'potential' only)
241 //! Distance between two angles in degrees (for fit type 'potential' only)
243 //! Slab distance (nm)
245 //! Minimum value the gaussian must have so that the force is actually evaluated
247 //! Additive constant for radial motion2 and flexible2 potentials (nm^2)
253 //! Number of rotation groups
255 //! Output frequency for main rotation outfile
257 //! Output frequency for per-slab data
265 //! Number of interactive atoms
267 //! The global indices of the interactive atoms
273 //! Name of the swap group, e.g. NA, CL, SOL
275 //! Number of atoms in this group
277 //! The global ion group atoms numbers
279 //! Requested number of molecules of this type per compartment
280 int nmolReq[eCompNR];
285 //! Period between when a swap is attempted
287 //! Use mass-weighted positions in split group
288 gmx_bool massw_split[2];
289 /*! \brief Split cylinders defined by radius, upper and lower
290 * extension. The split cylinders define the channels and are
291 * each anchored in the center of the split group */
297 //! Coupling constant (number of swap attempt steps)
299 //! Ion counts may deviate from the requested values by +-threshold before a swap is done
301 //! Offset of the swap layer (='bulk') with respect to the compartment-defining layers
302 real bulkOffset[eCompNR];
303 //! Number of groups to be controlled
305 //! All swap groups, including split and solvent
309 struct t_inputrec // NOLINT (clang-analyzer-optin.performance.Padding)
312 explicit t_inputrec(const t_inputrec&) = delete;
313 t_inputrec& operator=(const t_inputrec&) = delete;
316 //! Integration method
317 IntegrationAlgorithm eI;
318 //! Number of steps to be taken
320 //! Used in checkpointing to separate chunks
322 //! Start at a stepcount >0 (used w. convert-tpr)
324 //! Frequency of energy calc. and T/P coupl. upd.
326 //! Group or verlet cutoffs
327 CutoffScheme cutoff_scheme;
328 //! Number of steps before pairlist is generated
330 //! Number of steps after which center of mass motion is removed
332 //! Center of mass motion removal algorithm
333 ComRemovalAlgorithm comm_mode;
334 //! Number of steps after which print to logfile
336 //! Number of steps after which X is output
338 //! Number of steps after which V is output
340 //! Number of steps after which F is output
342 //! Number of steps after which energies printed
344 //! Number of steps after which compressed trj (.xtc,.tng) is output
345 int nstxout_compressed;
346 //! Initial time (ps)
350 //! Whether we use multiple time stepping
352 //! The multiple time stepping levels
353 std::vector<gmx::MtsLevel> mtsLevels;
354 //! Precision of x in compressed trajectory file
355 real x_compression_precision;
356 //! Requested fourier_spacing, when nk? not set
357 real fourier_spacing;
358 //! Number of k vectors in x dimension for fourier methods for long range electrost.
360 //! Number of k vectors in y dimension for fourier methods for long range electrost.
362 //! Number of k vectors in z dimension for fourier methods for long range electrost.
364 //! Interpolation order for PME
366 //! Real space tolerance for Ewald, determines the real/reciprocal space relative weight
368 //! Real space tolerance for LJ-Ewald
370 //! Normal/3D ewald, or pseudo-2D LR corrections
371 EwaldGeometry ewald_geometry;
372 //! Epsilon for PME dipole correction
373 real epsilon_surface;
374 //! Type of combination rule in LJ-PME
375 LongRangeVdW ljpme_combination_rule;
376 //! Type of periodic boundary conditions
378 //! Periodic molecules
380 //! Continuation run: starting state is correct (ie. constrained)
381 gmx_bool bContinuation;
382 //! Temperature coupling
383 TemperatureCoupling etc;
384 //! Interval in steps for temperature coupling
386 //! Whether to print nose-hoover chains
387 gmx_bool bPrintNHChains;
388 //! Pressure coupling
389 PressureCoupling epc;
390 //! Pressure coupling type
391 PressureCouplingType epct;
392 //! Interval in steps for pressure coupling
394 //! Pressure coupling time (ps)
396 //! Reference pressure (kJ/(mol nm^3))
398 //! Compressibility ((mol nm^3)/kJ)
400 //! How to scale absolute reference coordinates
401 RefCoordScaling refcoord_scaling;
402 //! The COM of the posres atoms
404 //! The B-state COM of the posres atoms
406 //! Random seed for Andersen thermostat (obsolete)
408 //! Per atom pair energy drift tolerance (kJ/mol/ps/atom) for list buffer
410 //! Short range pairlist cut-off (nm)
412 //! Radius for test particle insertion
414 //! Type of electrostatics treatment
415 CoulombInteractionType coulombtype;
416 //! Modify the Coulomb interaction
417 InteractionModifiers coulomb_modifier;
418 //! Coulomb switch range start (nm)
419 real rcoulomb_switch;
420 //! Coulomb cutoff (nm)
422 //! Relative dielectric constant
424 //! Relative dielectric constant of the RF
426 //! Always false (no longer supported)
427 bool implicit_solvent;
428 //! Type of Van der Waals treatment
429 VanDerWaalsType vdwtype;
430 //! Modify the Van der Waals interaction
431 InteractionModifiers vdw_modifier;
432 //! Van der Waals switch range start (nm)
434 //! Van der Waals cutoff (nm)
436 //! Perform Long range dispersion corrections
437 DispersionCorrectionType eDispCorr;
438 //! Extension of the table beyond the cut-off, as well as the table length for 1-4 interac.
440 //! Tolerance for shake
442 //! Free energy calculations
443 FreeEnergyPerturbationType efep;
444 //! Data for the FEP state
445 std::unique_ptr<t_lambda> fepvals;
446 //! Whether to do simulated tempering
448 //! Variables for simulated tempering
449 std::unique_ptr<t_simtemp> simtempvals;
450 //! Whether expanded ensembles are used
452 //! Expanded ensemble parameters
453 std::unique_ptr<t_expanded> expandedvals;
454 //! Type of distance restraining
455 DistanceRestraintRefinement eDisre;
456 //! Force constant for time averaged distance restraints
458 //! Type of weighting of pairs in one restraints
459 DistanceRestraintWeighting eDisreWeighting;
460 //! Use combination of time averaged and instantaneous violations
461 gmx_bool bDisreMixed;
462 //! Frequency of writing pair distances to enx
464 //! Time constant for memory function in disres
466 //! Force constant for orientational restraints
468 //! Time constant for memory function in orires
470 //! Frequency of writing tr(SD) to energy output
472 //! The stepsize for updating
476 //! Number of iterations for convergence of steepest descent in relax_shells
478 //! Stepsize for directional minimization in relax_shells
480 //! Number of steps after which a steepest descents step is done while doing cg
482 //! Number of corrections to the Hessian to keep
484 //! Type of constraint algorithm
485 ConstraintAlgorithm eConstrAlg;
486 //! Order of the LINCS Projection Algorithm
488 //! Warn if any bond rotates more than this many degrees
490 //! Number of iterations in the final LINCS step
492 //! Use successive overrelaxation for shake
494 //! Friction coefficient for BD (amu/ps)
496 //! Random seed for SD and BD
498 //! The number of walls
500 //! The type of walls
502 //! The potentail is linear for r<=wall_r_linpot
504 //! The atom type for walls
505 int wall_atomtype[2];
506 //! Number density for walls
507 real wall_density[2];
508 //! Scaling factor for the box for Ewald
509 real wall_ewald_zfac;
511 /* COM pulling data */
512 //! Do we do COM pulling?
514 //! The data for center of mass pulling
515 std::unique_ptr<pull_params_t> pull;
518 //! Whether to use AWH biasing for PMF calculations
520 //! AWH biasing parameters
521 std::unique_ptr<gmx::AwhParams> awhParams;
523 /* Enforced rotation data */
524 //! Whether to calculate enforced rotation potential(s)
526 //! The data for enforced rotation potentials
529 //! Whether to do ion/water position exchanges (CompEL)
530 SwapType eSwapCoords;
531 //! Swap data structure.
534 //! Whether the tpr makes an interactive MD session possible.
536 //! Interactive molecular dynamics
539 //! Acceleration for viscosity calculation
541 //! Triclinic deformation velocities (nm/ps)
543 /*! \brief User determined parameters */
556 //! QM/MM calculation
559 /* Fields for removed features go here (better caching) */
560 //! Whether AdResS is enabled - always false if a valid .tpr was read
562 //! Whether twin-range scheme is active - always false if a valid .tpr was read
563 gmx_bool useTwinRange;
564 //! Whether we have constant acceleration - removed in GROMACS 2022
565 bool useConstantAcceleration;
567 //! KVT object that contains input parameters converted to the new style.
568 gmx::KeyValueTreeObject* params;
570 //! KVT for storing simulation parameters that are not part of the mdp file.
571 std::unique_ptr<gmx::KeyValueTreeObject> internalParameters;
574 int ir_optimal_nstcalcenergy(const t_inputrec* ir);
576 int tcouple_min_integration_steps(TemperatureCoupling etc);
578 int ir_optimal_nsttcouple(const t_inputrec* ir);
580 int pcouple_min_integration_steps(PressureCoupling epc);
582 int ir_optimal_nstpcouple(const t_inputrec* ir);
584 /* Returns if the Coulomb force or potential is switched to zero */
585 gmx_bool ir_coulomb_switched(const t_inputrec* ir);
587 /* Returns if the Coulomb interactions are zero beyond the rcoulomb.
588 * Note: always returns TRUE for the Verlet cut-off scheme.
590 gmx_bool ir_coulomb_is_zero_at_cutoff(const t_inputrec* ir);
592 /* As ir_coulomb_is_zero_at_cutoff, but also returns TRUE for user tabulated
593 * interactions, since these might be zero beyond rcoulomb.
595 gmx_bool ir_coulomb_might_be_zero_at_cutoff(const t_inputrec* ir);
597 /* Returns if the Van der Waals force or potential is switched to zero */
598 gmx_bool ir_vdw_switched(const t_inputrec* ir);
600 /* Returns if the Van der Waals interactions are zero beyond the rvdw.
601 * Note: always returns TRUE for the Verlet cut-off scheme.
603 gmx_bool ir_vdw_is_zero_at_cutoff(const t_inputrec* ir);
605 /* As ir_vdw_is_zero_at_cutoff, but also returns TRUE for user tabulated
606 * interactions, since these might be zero beyond rvdw.
608 gmx_bool ir_vdw_might_be_zero_at_cutoff(const t_inputrec* ir);
610 /*! \brief Free memory from input record.
612 * All arrays and pointers will be freed.
614 * \param[in] ir The data structure
616 void done_inputrec(t_inputrec* ir);
618 void pr_inputrec(FILE* fp, int indent, const char* title, const t_inputrec* ir, gmx_bool bMDPformat);
620 void cmp_inputrec(FILE* fp, const t_inputrec* ir1, const t_inputrec* ir2, real ftol, real abstol);
622 void comp_pull_AB(FILE* fp, const pull_params_t& pull, real ftol, real abstol);
625 gmx_bool inputrecDeform(const t_inputrec* ir);
627 gmx_bool inputrecDynamicBox(const t_inputrec* ir);
629 gmx_bool inputrecPreserveShape(const t_inputrec* ir);
631 gmx_bool inputrecNeedMutot(const t_inputrec* ir);
633 gmx_bool inputrecTwinRange(const t_inputrec* ir);
635 gmx_bool inputrecExclForces(const t_inputrec* ir);
637 gmx_bool inputrecNptTrotter(const t_inputrec* ir);
639 gmx_bool inputrecNvtTrotter(const t_inputrec* ir);
641 gmx_bool inputrecNphTrotter(const t_inputrec* ir);
643 /*! \brief Return true if the simulation is 2D periodic with two walls. */
644 bool inputrecPbcXY2Walls(const t_inputrec* ir);
646 /*! \brief Returns true for MD integator with T and/or P-coupling that supports
647 * calculating a conserved energy quantity.
649 * Note that dynamical integrators without T and P coupling (ie NVE)
650 * return false, i.e. the return value refers to whether there
651 * is a conserved quantity different than the total energy.
653 bool integratorHasConservedEnergyQuantity(const t_inputrec* ir);
655 /*! \brief Returns true when temperature is coupled or constant. */
656 bool integratorHasReferenceTemperature(const t_inputrec* ir);
658 /*! \brief Return the number of bounded dimensions
660 * \param[in] ir The input record with MD parameters
661 * \return the number of dimensions in which
662 * the coordinates of the particles are bounded, starting at X.
664 int inputrec2nboundeddim(const t_inputrec* ir);
666 /*! \brief Returns the number of degrees of freedom in center of mass motion
668 * \param[in] ir The inputrec structure
669 * \return the number of degrees of freedom of the center of mass
671 int ndof_com(const t_inputrec* ir);
673 /*! \brief Returns the maximum reference temperature over all coupled groups
675 * Returns 0 for energy minimization and normal mode computation.
676 * Returns -1 for MD without temperature coupling.
678 * \param[in] ir The inputrec structure
680 real maxReferenceTemperature(const t_inputrec& ir);
682 /*! \brief Returns whether there is an Ewald surface contribution
684 bool haveEwaldSurfaceContribution(const t_inputrec& ir);
686 /*! \brief Check if calculation of the specific FEP type was requested.
688 * \param[in] ir Input record.
689 * \param[in] fepType Free-energy perturbation type to check for.
691 * \returns If the \p fepType is perturbed in this run.
693 bool haveFreeEnergyType(const t_inputrec& ir, int fepType);
695 #endif /* GMX_MDTYPES_INPUTREC_H */