*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014,2015,2016,2017, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team.
+ * Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
#ifndef GMX_MDTYPES_INPUTREC_H
#define GMX_MDTYPES_INPUTREC_H
-#include <stdio.h>
+#include <cstdio>
+
+#include <memory>
+#include <vector>
#include "gromacs/math/vectypes.h"
#include "gromacs/mdtypes/md_enums.h"
-#include "gromacs/utility/basedefinitions.h"
+#include "gromacs/utility/enumerationhelpers.h"
#include "gromacs/utility/real.h"
-#define EGP_EXCL (1<<0)
-#define EGP_TABLE (1<<1)
+#define EGP_EXCL (1 << 0)
+#define EGP_TABLE (1 << 1)
+struct gmx_enfrot;
+struct gmx_enfrotgrp;
struct pull_params_t;
namespace gmx
{
+class Awh;
+class AwhParams;
class KeyValueTreeObject;
-}
-
-typedef struct t_grpopts {
- int ngtc; /* # T-Coupl groups */
- int nhchainlength; /* # of nose-hoover chains per group */
- int ngacc; /* # Accelerate groups */
- int ngfrz; /* # Freeze groups */
- int ngener; /* # Ener groups */
- real *nrdf; /* Nr of degrees of freedom in a group */
- real *ref_t; /* Coupling temperature per group */
- int *annealing; /* No/simple/periodic SA for each group */
- int *anneal_npoints; /* Number of annealing time points per grp */
- real **anneal_time; /* For ea. group: Time points */
- real **anneal_temp; /* For ea. grp: Temperature at these times */
- /* Final temp after all intervals is ref_t */
- real *tau_t; /* Tau coupling time */
- rvec *acc; /* Acceleration per group */
- ivec *nFreeze; /* Freeze the group in each direction ? */
- int *egp_flags; /* Exclusions/tables of energy group pairs */
+struct MtsLevel;
+} // namespace gmx
+
+struct t_grpopts
+{
+ //! Number of T-Coupl groups
+ int ngtc = 0;
+ //! Number of of Nose-Hoover chains per group
+ int nhchainlength = 0;
+ //! Number of Accelerate groups
+ int ngacc;
+ //! Number of Freeze groups
+ int ngfrz = 0;
+ //! Number of Energy groups
+ int ngener = 0;
+ //! Number of degrees of freedom in a group
+ real* nrdf = nullptr;
+ //! Coupling temperature per group
+ real* ref_t = nullptr;
+ //! No/simple/periodic simulated annealing for each group
+ SimulatedAnnealing* annealing = nullptr;
+ //! Number of annealing time points per group
+ int* anneal_npoints = nullptr;
+ //! For each group: Time points
+ real** anneal_time = nullptr;
+ //! For each group: Temperature at these times. Final temp after all intervals is ref_t
+ real** anneal_temp = nullptr;
+ //! Tau coupling time
+ real* tau_t = nullptr;
+ //! Acceleration per group
+ rvec* acceleration = nullptr;
+ //! Whether the group will be frozen in each direction
+ ivec* nFreeze = nullptr;
+ //! Exclusions/tables of energy group pairs
+ int* egp_flags = nullptr;
/* QMMM stuff */
- int ngQM; /* nr of QM groups */
- int *QMmethod; /* Level of theory in the QM calculation */
- int *QMbasis; /* Basisset in the QM calculation */
- int *QMcharge; /* Total charge in the QM region */
- int *QMmult; /* Spin multiplicicty in the QM region */
- gmx_bool *bSH; /* surface hopping (diabatic hop only) */
- int *CASorbitals; /* number of orbiatls in the active space */
- int *CASelectrons; /* number of electrons in the active space */
- real *SAon; /* at which gap (A.U.) the SA is switched on */
- real *SAoff;
- int *SAsteps; /* in how many steps SA goes from 1-1 to 0.5-0.5*/
-} t_grpopts;
-
-typedef struct t_simtemp {
- int eSimTempScale; /* simulated temperature scaling; linear or exponential */
- real simtemp_low; /* the low temperature for simulated tempering */
- real simtemp_high; /* the high temperature for simulated tempering */
- real *temperatures; /* the range of temperatures used for simulated tempering */
-} t_simtemp;
-
-typedef struct t_lambda {
- int nstdhdl; /* The frequency for calculating dhdl */
- double init_lambda; /* 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) */
- int init_fep_state; /* the initial number of the state */
- double delta_lambda; /* change of lambda per time step (fraction of (0.1) */
- int edHdLPrintEnergy; /* print no, total or potential energies in dhdl */
- int n_lambda; /* The number of foreign lambda points */
- double **all_lambda; /* The array of all lambda values */
- int lambda_neighbors; /* The number of neighboring lambda states to
- calculate the energy for in up and down directions
- (-1 for all) */
- int lambda_start_n; /* The first lambda to calculate energies for */
- int lambda_stop_n; /* The last lambda +1 to calculate energies for */
- real sc_alpha; /* free energy soft-core parameter */
- int sc_power; /* lambda power for soft-core interactions */
- real sc_r_power; /* r power for soft-core interactions */
- real sc_sigma; /* free energy soft-core sigma when c6 or c12=0 */
- real sc_sigma_min; /* free energy soft-core sigma for ????? */
- gmx_bool bScCoul; /* use softcore for the coulomb portion as well (default FALSE) */
- gmx_bool separate_dvdl[efptNR]; /* whether to print the dvdl term associated with
- this term; if it is not specified as separate,
- it is lumped with the FEP term */
- int separate_dhdl_file; /* whether to write a separate dhdl.xvg file
- note: NOT a gmx_bool, but an enum */
- int dhdl_derivatives; /* whether to calculate+write dhdl derivatives
- note: NOT a gmx_bool, but an enum */
- int dh_hist_size; /* The maximum table size for the dH histogram */
- double dh_hist_spacing; /* The spacing for the dH histogram */
-} t_lambda;
-
-typedef struct t_expanded {
- int nstexpanded; /* The frequency of expanded ensemble state changes */
- int elamstats; /* which type of move updating do we use for lambda monte carlo (or no for none) */
- int elmcmove; /* what move set will be we using for state space moves */
- int elmceq; /* the method we use to decide of we have equilibrated the weights */
- int equil_n_at_lam; /* the minumum number of samples at each lambda for deciding whether we have reached a minimum */
- real equil_wl_delta; /* WL delta at which we stop equilibrating weights */
- real equil_ratio; /* use the ratio of weights (ratio of minimum to maximum) to decide when to stop equilibrating */
- int equil_steps; /* after equil_steps steps we stop equilibrating the weights */
- int equil_samples; /* after equil_samples total samples (steps/nstfep), we stop equilibrating the weights */
- int lmc_seed; /* random number seed for lambda mc switches */
- gmx_bool minvar; /* whether to use minumum variance weighting */
- int minvarmin; /* the number of samples needed before kicking into minvar routine */
- real minvar_const; /* the offset for the variance in MinVar */
- int c_range; /* range of cvalues used for BAR */
- gmx_bool bSymmetrizedTMatrix; /* whether to print symmetrized matrices */
- int nstTij; /* How frequently to print the transition matrices */
- int lmc_repeats; /* number of repetitions in the MC lambda jumps */ /*MRS -- VERIFY THIS */
- int lmc_forced_nstart; /* minimum number of samples for each state before free sampling */ /* MRS -- VERIFY THIS! */
- int gibbsdeltalam; /* distance in lambda space for the gibbs interval */
- real wl_scale; /* scaling factor for wang-landau */
- real wl_ratio; /* ratio between largest and smallest number for freezing the weights */
- real init_wl_delta; /* starting delta for wang-landau */
- gmx_bool bWLoneovert; /* use one over t convergece for wang-landau when the delta get sufficiently small */
- gmx_bool bInit_weights; /* did we initialize the weights? TODO: REMOVE FOR 5.0, no longer needed with new logic */
- real mc_temp; /* To override the main temperature, or define it if it's not defined */
- real *init_lambda_weights; /* user-specified initial weights to start with */
-} t_expanded;
-
-
-/* Abstract types for enforced rotation only defined in pull_rotation.c */
-typedef struct gmx_enfrot *gmx_enfrot_t;
-typedef struct gmx_enfrotgrp *gmx_enfrotgrp_t;
-
-typedef struct {
- int eType; /* Rotation type for this group */
- int bMassW; /* Use mass-weighed positions? */
- int nat; /* Number of atoms in the group */
- int *ind; /* The global atoms numbers */
- rvec *x_ref; /* The reference positions */
- rvec vec; /* The normalized rotation vector */
- real rate; /* Rate of rotation (degree/ps) */
- real k; /* Force constant (kJ/(mol nm^2) */
- rvec pivot; /* Pivot point of rotation axis (nm) */
- int eFittype; /* Type of fit to determine actual group angle */
- int PotAngle_nstep; /* Number of angles around the reference angle
- for which the rotation potential is also
- evaluated (for fit type 'potential' only) */
- real PotAngle_step; /* Distance between two angles in degrees (for
- fit type 'potential' only) */
- real slab_dist; /* Slab distance (nm) */
- real min_gaussian; /* Minimum value the gaussian must have so that
- the force is actually evaluated */
- real eps; /* Additive constant for radial motion2 and
- flexible2 potentials (nm^2) */
- gmx_enfrotgrp_t enfrotgrp; /* Stores non-inputrec rotation data per group */
-} t_rotgrp;
-
-typedef struct t_rot {
- int ngrp; /* Number of rotation groups */
- int nstrout; /* Output frequency for main rotation outfile */
- int nstsout; /* Output frequency for per-slab data */
- t_rotgrp *grp; /* Groups to rotate */
- gmx_enfrot_t enfrot; /* Stores non-inputrec enforced rotation data */
-} t_rot;
-
-/* Abstract type for IMD only defined in IMD.c */
-struct t_gmx_IMD;
-
-typedef struct t_IMD {
- int nat; /* Number of interactive atoms */
- int *ind; /* The global indices of the interactive atoms */
- struct t_gmx_IMD *setup; /* Stores non-inputrec IMD data */
-} t_IMD;
-
-/* Abstract types for position swapping only defined in swapcoords.cpp */
-typedef struct t_swap *gmx_swapcoords_t;
-
-typedef struct t_swapGroup {
- char *molname; /* Name of the swap group, e.g. NA, CL, SOL */
- int nat; /* Number of atoms in this group */
- int *ind; /* The global ion group atoms numbers */
- int nmolReq[eCompNR]; /* Requested number of molecules of this type
- per compartment */
-} t_swapGroup;
-
-typedef struct t_swapcoords {
- int nstswap; /* Every how many steps a swap is attempted? */
- gmx_bool massw_split[2]; /* Use mass-weighted positions in split group? */
- real cyl0r, cyl1r; /* Split cylinders defined by radius, upper and */
- real cyl0u, cyl1u; /* ... lower extension. The split cylinders de- */
- real cyl0l, cyl1l; /* ... fine the channels and are each anchored */
- /* ... in the center of the split group */
- int nAverage; /* Coupling constant (nr of swap attempt steps) */
- real threshold; /* Ion counts may deviate from the requested
- values by +-threshold before a swap is done */
- real bulkOffset[eCompNR]; /* Offset of the swap layer (='bulk') w.r.t.
- the compartment-defining layers */
- int ngrp; /* Number of groups to be controlled */
- t_swapGroup *grp; /* All swap groups, including split and solvent */
- gmx_swapcoords_t si_priv; /* swap private data accessible in
- * swapcoords.cpp */
-} t_swapcoords;
-
-struct t_inputrec
+ //! Number of QM groups
+ int ngQM = 0;
+};
+
+struct t_simtemp
+{
+ //! Simulated temperature scaling; linear or exponential
+ SimulatedTempering eSimTempScale = SimulatedTempering::Default;
+ //! The low temperature for simulated tempering
+ real simtemp_low = 0;
+ //! The high temperature for simulated tempering
+ real simtemp_high = 0;
+ //! The range of temperatures used for simulated tempering
+ std::vector<real> temperatures;
+};
+
+struct t_lambda
+{
+ //! The frequency for calculating dhdl
+ int nstdhdl = 0;
+ //! 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)
+ double init_lambda = 0;
+ //! The initial number of the state
+ int init_fep_state = 0;
+ //! Change of lambda per time step (fraction of (0.1)
+ double delta_lambda = 0;
+ //! Print no, total or potential energies in dhdl
+ FreeEnergyPrintEnergy edHdLPrintEnergy = FreeEnergyPrintEnergy::Default;
+ //! The number of foreign lambda points
+ int n_lambda = 0;
+ //! The array of all lambda values
+ gmx::EnumerationArray<FreeEnergyPerturbationCouplingType, std::vector<double>> all_lambda;
+ //! The number of neighboring lambda states to calculate the energy for in up and down directions (-1 for all)
+ int lambda_neighbors = 0;
+ //! The first lambda to calculate energies for
+ int lambda_start_n = 0;
+ //! The last lambda +1 to calculate energies for
+ int lambda_stop_n = 0;
+ //! Free energy soft-core parameter
+ real sc_alpha = 0;
+ //! Lambda power for soft-core interactions
+ int sc_power = 0;
+ //! R power for soft-core interactions
+ real sc_r_power = 0;
+ //! Free energy soft-core sigma when c6 or c12=0
+ real sc_sigma = 0;
+ //! Free energy soft-core sigma for ?????
+ real sc_sigma_min = 0;
+ //! Use softcore for the coulomb portion as well (default FALSE)
+ bool bScCoul = false;
+ //! The specific soft-core function to use
+ SoftcoreType softcoreFunction = SoftcoreType::Beutler;
+ //! scale for the linearization point for the vdw interaction with gapsys soft-core
+ real scGapsysScaleLinpointLJ = 0.85;
+ //! scale for the linearization point for the coulomb interaction with gapsys soft-core
+ real scGapsysScaleLinpointQ = 0.3;
+ //! lower bound for c12/c6 in gapsys soft-core
+ real scGapsysSigmaLJ = 0.3;
+ //! Whether to print the dvdl term associated with this term; if it is not specified as separate, it is lumped with the FEP term
+ gmx::EnumerationArray<FreeEnergyPerturbationCouplingType, bool> separate_dvdl;
+ //! Whether to write a separate dhdl.xvg file note: NOT a gmx_bool, but an enum
+ SeparateDhdlFile separate_dhdl_file = SeparateDhdlFile::Default;
+ //! Whether to calculate+write dhdl derivatives note: NOT a gmx_bool, but an enum
+ DhDlDerivativeCalculation dhdl_derivatives = DhDlDerivativeCalculation::Default;
+ //! The maximum table size for the dH histogram
+ int dh_hist_size = 0;
+ //! The spacing for the dH histogram
+ double dh_hist_spacing = 0;
+};
+
+struct t_expanded
+{
+ //! The frequency of expanded ensemble state changes
+ int nstexpanded = 0;
+ //! Which type of move updating do we use for lambda monte carlo (or no for none)
+ LambdaWeightCalculation elamstats = LambdaWeightCalculation::Default;
+ //! What move set will be we using for state space moves
+ LambdaMoveCalculation elmcmove = LambdaMoveCalculation::Default;
+ //! The method we use to decide of we have equilibrated the weights
+ LambdaWeightWillReachEquilibrium elmceq = LambdaWeightWillReachEquilibrium::Default;
+ //! The minumum number of samples at each lambda for deciding whether we have reached a minimum
+ int equil_n_at_lam = 0;
+ //! Wang-Landau delta at which we stop equilibrating weights
+ real equil_wl_delta = 0;
+ //! Use the ratio of weights (ratio of minimum to maximum) to decide when to stop equilibrating
+ real equil_ratio = 0;
+ //! After equil_steps steps we stop equilibrating the weights
+ int equil_steps = 0;
+ //! After equil_samples total samples (steps/nstfep), we stop equilibrating the weights
+ int equil_samples = 0;
+ //! Random number seed for lambda mc switches
+ int lmc_seed = 0;
+ //! Whether to use minumum variance weighting
+ bool minvar = false;
+ //! The number of samples needed before kicking into minvar routine
+ int minvarmin = 0;
+ //! The offset for the variance in MinVar
+ real minvar_const = 0;
+ //! Range of cvalues used for BAR
+ int c_range = 0;
+ //! Whether to print symmetrized matrices
+ bool bSymmetrizedTMatrix = false;
+ //! How frequently to print the transition matrices
+ int nstTij = 0;
+ //! Number of repetitions in the MC lambda jumps MRS -- VERIFY THIS
+ int lmc_repeats = 0;
+ //! Minimum number of samples for each state before free sampling MRS -- VERIFY THIS!
+ int lmc_forced_nstart = 0;
+ //! Distance in lambda space for the gibbs interval
+ int gibbsdeltalam = 0;
+ //! Scaling factor for Wang-Landau
+ real wl_scale = 0;
+ //! Ratio between largest and smallest number for freezing the weights
+ real wl_ratio = 0;
+ //! Starting delta for Wang-Landau
+ real init_wl_delta = 0;
+ //! Use one over t convergence for Wang-Landau when the delta get sufficiently small
+ bool bWLoneovert = false;
+ //! Did we initialize the weights? TODO: REMOVE FOR 5.0, no longer needed with new logic
+ bool bInit_weights = false;
+ //! To override the main temperature, or define it if it's not defined
+ real mc_temp = 0;
+ //! User-specified initial weights to start with
+ std::vector<real> init_lambda_weights;
+};
+
+struct t_rotgrp
+{
+ //! Rotation type for this group
+ EnforcedRotationGroupType eType;
+ //! Use mass-weighed positions?
+ bool bMassW;
+ //! Number of atoms in the group
+ int nat;
+ //! The global atoms numbers
+ int* ind;
+ //! The reference positions
+ rvec* x_ref;
+ //! The normalized rotation vector
+ rvec inputVec;
+ //! Rate of rotation (degree/ps)
+ real rate;
+ //! Force constant (kJ/(mol nm^2)
+ real k;
+ //! Pivot point of rotation axis (nm)
+ rvec pivot;
+ //! Type of fit to determine actual group angle
+ RotationGroupFitting eFittype;
+ //! Number of angles around the reference angle for which the rotation potential is also evaluated (for fit type 'potential' only)
+ int PotAngle_nstep;
+ //! Distance between two angles in degrees (for fit type 'potential' only)
+ real PotAngle_step;
+ //! Slab distance (nm)
+ real slab_dist;
+ //! Minimum value the gaussian must have so that the force is actually evaluated
+ real min_gaussian;
+ //! Additive constant for radial motion2 and flexible2 potentials (nm^2)
+ real eps;
+};
+
+struct t_rot
+{
+ //! Number of rotation groups
+ int ngrp;
+ //! Output frequency for main rotation outfile
+ int nstrout;
+ //! Output frequency for per-slab data
+ int nstsout;
+ //! Groups to rotate
+ t_rotgrp* grp;
+};
+
+struct t_IMD
+{
+ //! Number of interactive atoms
+ int nat;
+ //! The global indices of the interactive atoms
+ int* ind;
+};
+
+struct t_swapGroup
+{
+ //! Name of the swap group, e.g. NA, CL, SOL
+ char* molname;
+ //! Number of atoms in this group
+ int nat;
+ //! The global ion group atoms numbers
+ int* ind;
+ //! Requested number of molecules of this type per compartment
+ gmx::EnumerationArray<Compartment, int> nmolReq;
+};
+
+struct t_swapcoords
+{
+ //! Period between when a swap is attempted
+ int nstswap;
+ //! Use mass-weighted positions in split group
+ bool massw_split[2];
+ /*! \brief Split cylinders defined by radius, upper and lower
+ * extension. The split cylinders define the channels and are
+ * each anchored in the center of the split group */
+ /**@{*/
+ real cyl0r, cyl1r;
+ real cyl0u, cyl1u;
+ real cyl0l, cyl1l;
+ /**@}*/
+ //! Coupling constant (number of swap attempt steps)
+ int nAverage;
+ //! Ion counts may deviate from the requested values by +-threshold before a swap is done
+ real threshold;
+ //! Offset of the swap layer (='bulk') with respect to the compartment-defining layers
+ gmx::EnumerationArray<Compartment, real> bulkOffset;
+ //! Number of groups to be controlled
+ int ngrp;
+ //! All swap groups, including split and solvent
+ t_swapGroup* grp;
+};
+
+struct t_inputrec // NOLINT (clang-analyzer-optin.performance.Padding)
{
t_inputrec();
- explicit t_inputrec(const t_inputrec &) = delete;
- t_inputrec &operator=(const t_inputrec &) = delete;
+ explicit t_inputrec(const t_inputrec&) = delete;
+ t_inputrec& operator=(const t_inputrec&) = delete;
~t_inputrec();
- int eI; /* Integration method */
- gmx_int64_t nsteps; /* number of steps to be taken */
- int simulation_part; /* Used in checkpointing to separate chunks */
- gmx_int64_t init_step; /* start at a stepcount >0 (used w. convert-tpr) */
- int nstcalcenergy; /* frequency of energy calc. and T/P coupl. upd. */
- int cutoff_scheme; /* group or verlet cutoffs */
- int ns_type; /* which ns method should we use? */
- int nstlist; /* number of steps before pairlist is generated */
- int ndelta; /* number of cells per rlong */
- int nstcomm; /* number of steps after which center of mass */
- /* motion is removed */
- int comm_mode; /* Center of mass motion removal algorithm */
- int nstlog; /* number of steps after which print to logfile */
- int nstxout; /* number of steps after which X is output */
- int nstvout; /* id. for V */
- int nstfout; /* id. for F */
- int nstenergy; /* number of steps after which energies printed */
- int nstxout_compressed; /* id. for compressed trj (.xtc,.tng) */
- double init_t; /* initial time (ps) */
- double delta_t; /* time step (ps) */
- real x_compression_precision; /* precision of x in compressed trajectory file */
- real fourier_spacing; /* requested fourier_spacing, when nk? not set */
- int nkx, nky, nkz; /* number of k vectors in each spatial dimension*/
- /* for fourier methods for long range electrost.*/
- int pme_order; /* interpolation order for PME */
- real ewald_rtol; /* Real space tolerance for Ewald, determines */
- /* the real/reciprocal space relative weight */
- real ewald_rtol_lj; /* Real space tolerance for LJ-Ewald */
- int ewald_geometry; /* normal/3d ewald, or pseudo-2d LR corrections */
- real epsilon_surface; /* Epsilon for PME dipole correction */
- int ljpme_combination_rule; /* Type of combination rule in LJ-PME */
- int ePBC; /* Type of periodic boundary conditions */
- int bPeriodicMols; /* Periodic molecules */
- gmx_bool bContinuation; /* Continuation run: starting state is correct */
- int etc; /* temperature coupling */
- int nsttcouple; /* interval in steps for temperature coupling */
- gmx_bool bPrintNHChains; /* whether to print nose-hoover chains */
- int epc; /* pressure coupling */
- int epct; /* pressure coupling type */
- int nstpcouple; /* interval in steps for pressure coupling */
- real tau_p; /* pressure coupling time (ps) */
- tensor ref_p; /* reference pressure (kJ/(mol nm^3)) */
- tensor compress; /* compressability ((mol nm^3)/kJ) */
- int refcoord_scaling; /* How to scale absolute reference coordinates */
- rvec posres_com; /* The COM of the posres atoms */
- rvec posres_comB; /* The B-state COM of the posres atoms */
- int andersen_seed; /* Random seed for Andersen thermostat (obsolete) */
- real verletbuf_tol; /* Per atom pair energy drift tolerance (kJ/mol/ps/atom) for list buffer */
- real rlist; /* short range pairlist cut-off (nm) */
- real rtpi; /* Radius for test particle insertion */
- int coulombtype; /* Type of electrostatics treatment */
- int coulomb_modifier; /* Modify the Coulomb interaction */
- real rcoulomb_switch; /* Coulomb switch range start (nm) */
- real rcoulomb; /* Coulomb cutoff (nm) */
- real epsilon_r; /* relative dielectric constant */
- real epsilon_rf; /* relative dielectric constant of the RF */
- int implicit_solvent; /* No (=explicit water), or GBSA solvent models */
- int gb_algorithm; /* Algorithm to use for calculation Born radii */
- int nstgbradii; /* Frequency of updating Generalized Born radii */
- real rgbradii; /* Cutoff for GB radii calculation */
- real gb_saltconc; /* Salt concentration (M) for GBSA models */
- real gb_epsilon_solvent; /* dielectric coeff. of implicit solvent */
- real gb_obc_alpha; /* 1st scaling factor for Bashford-Case GB */
- real gb_obc_beta; /* 2nd scaling factor for Bashford-Case GB */
- real gb_obc_gamma; /* 3rd scaling factor for Bashford-Case GB */
- real gb_dielectric_offset; /* Dielectric offset for Still/HCT/OBC */
- int sa_algorithm; /* Algorithm for SA part of GBSA */
- real sa_surface_tension; /* Energy factor for SA part of GBSA */
- int vdwtype; /* Type of Van der Waals treatment */
- int vdw_modifier; /* Modify the VdW interaction */
- real rvdw_switch; /* Van der Waals switch range start (nm) */
- real rvdw; /* Van der Waals cutoff (nm) */
- int eDispCorr; /* Perform Long range dispersion corrections */
- real tabext; /* Extension of the table beyond the cut-off, *
- * as well as the table length for 1-4 interac. */
- real shake_tol; /* tolerance for shake */
- int efep; /* free energy calculations */
- t_lambda *fepvals; /* Data for the FEP state */
- gmx_bool bSimTemp; /* Whether to do simulated tempering */
- t_simtemp *simtempvals; /* Variables for simulated tempering */
- gmx_bool bExpanded; /* Whether expanded ensembles are used */
- t_expanded *expandedvals; /* Expanded ensemble parameters */
- int eDisre; /* Type of distance restraining */
- real dr_fc; /* force constant for ta_disre */
- int eDisreWeighting; /* type of weighting of pairs in one restraints */
- gmx_bool bDisreMixed; /* Use comb of time averaged and instan. viol's */
- int nstdisreout; /* frequency of writing pair distances to enx */
- real dr_tau; /* time constant for memory function in disres */
- real orires_fc; /* force constant for orientational restraints */
- real orires_tau; /* time constant for memory function in orires */
- int nstorireout; /* frequency of writing tr(SD) to enx */
- real em_stepsize; /* The stepsize for updating */
- real em_tol; /* The tolerance */
- int niter; /* Number of iterations for convergence of */
- /* steepest descent in relax_shells */
- real fc_stepsize; /* Stepsize for directional minimization */
- /* in relax_shells */
- int nstcgsteep; /* number of steps after which a steepest */
- /* descents step is done while doing cg */
- int nbfgscorr; /* Number of corrections to the hessian to keep */
- int eConstrAlg; /* Type of constraint algorithm */
- int nProjOrder; /* Order of the LINCS Projection Algorithm */
- real LincsWarnAngle; /* If bond rotates more than %g degrees, warn */
- int nLincsIter; /* Number of iterations in the final Lincs step */
- gmx_bool bShakeSOR; /* Use successive overrelaxation for shake */
- real bd_fric; /* Friction coefficient for BD (amu/ps) */
- gmx_int64_t ld_seed; /* Random seed for SD and BD */
- int nwall; /* The number of walls */
- int wall_type; /* The type of walls */
- real wall_r_linpot; /* The potentail is linear for r<=wall_r_linpot */
- int wall_atomtype[2]; /* The atom type for walls */
- real wall_density[2]; /* Number density for walls */
- real wall_ewald_zfac; /* Scaling factor for the box for Ewald */
+ //! Integration method
+ IntegrationAlgorithm eI = IntegrationAlgorithm::Default;
+ //! Number of steps to be taken
+ int64_t nsteps = 0;
+ //! Used in checkpointing to separate chunks
+ int simulation_part = 0;
+ //! Start at a stepcount >0 (used w. convert-tpr)
+ int64_t init_step = 0;
+ //! Frequency of energy calc. and T/P coupl. upd.
+ int nstcalcenergy = 0;
+ //! Group or verlet cutoffs
+ CutoffScheme cutoff_scheme = CutoffScheme::Default;
+ //! Number of steps before pairlist is generated
+ int nstlist = 0;
+ //! Number of steps after which center of mass motion is removed
+ int nstcomm = 0;
+ //! Center of mass motion removal algorithm
+ ComRemovalAlgorithm comm_mode = ComRemovalAlgorithm::Default;
+ //! Number of steps after which print to logfile
+ int nstlog = 0;
+ //! Number of steps after which X is output
+ int nstxout = 0;
+ //! Number of steps after which V is output
+ int nstvout = 0;
+ //! Number of steps after which F is output
+ int nstfout = 0;
+ //! Number of steps after which energies printed
+ int nstenergy = 0;
+ //! Number of steps after which compressed trj (.xtc,.tng) is output
+ int nstxout_compressed = 0;
+ //! Initial time (ps)
+ double init_t = 0;
+ //! Time step (ps)
+ double delta_t = 0;
+ //! Whether we use multiple time stepping
+ bool useMts = false;
+ //! The multiple time stepping levels
+ std::vector<gmx::MtsLevel> mtsLevels;
+ //! Precision of x in compressed trajectory file
+ real x_compression_precision = 0;
+ //! Requested fourier_spacing, when nk? not set
+ real fourier_spacing = 0;
+ //! Number of k vectors in x dimension for fourier methods for long range electrost.
+ int nkx = 0;
+ //! Number of k vectors in y dimension for fourier methods for long range electrost.
+ int nky = 0;
+ //! Number of k vectors in z dimension for fourier methods for long range electrost.
+ int nkz = 0;
+ //! Interpolation order for PME
+ int pme_order = 0;
+ //! Real space tolerance for Ewald, determines the real/reciprocal space relative weight
+ real ewald_rtol = 0;
+ //! Real space tolerance for LJ-Ewald
+ real ewald_rtol_lj = 0;
+ //! Normal/3D ewald, or pseudo-2D LR corrections
+ EwaldGeometry ewald_geometry = EwaldGeometry::Default;
+ //! Epsilon for PME dipole correction
+ real epsilon_surface = 0;
+ //! Type of combination rule in LJ-PME
+ LongRangeVdW ljpme_combination_rule = LongRangeVdW::Default;
+ //! Type of periodic boundary conditions
+ PbcType pbcType = PbcType::Default;
+ //! Periodic molecules
+ bool bPeriodicMols = false;
+ //! Continuation run: starting state is correct (ie. constrained)
+ bool bContinuation = false;
+ //! Temperature coupling
+ TemperatureCoupling etc = TemperatureCoupling::Default;
+ //! Interval in steps for temperature coupling
+ int nsttcouple = 0;
+ //! Whether to print nose-hoover chains
+ bool bPrintNHChains = false;
+ //! Pressure coupling
+ PressureCoupling epc = PressureCoupling::Default;
+ //! Pressure coupling type
+ PressureCouplingType epct = PressureCouplingType::Default;
+ //! Interval in steps for pressure coupling
+ int nstpcouple = 0;
+ //! Pressure coupling time (ps)
+ real tau_p = 0;
+ //! Reference pressure (kJ/(mol nm^3))
+ tensor ref_p = { { 0 } };
+ //! Compressibility ((mol nm^3)/kJ)
+ tensor compress = { { 0 } };
+ //! How to scale absolute reference coordinates
+ RefCoordScaling refcoord_scaling = RefCoordScaling::Default;
+ //! The COM of the posres atoms
+ rvec posres_com = { 0, 0, 0 };
+ //! The B-state COM of the posres atoms
+ rvec posres_comB = { 0, 0, 0 };
+ //! Random seed for Andersen thermostat (obsolete)
+ int andersen_seed = 0;
+ //! Per atom pair energy drift tolerance (kJ/mol/ps/atom) for list buffer
+ real verletbuf_tol = 0;
+ //! Short range pairlist cut-off (nm)
+ real rlist = 0;
+ //! Radius for test particle insertion
+ real rtpi = 0;
+ //! Type of electrostatics treatment
+ CoulombInteractionType coulombtype = CoulombInteractionType::Default;
+ //! Modify the Coulomb interaction
+ InteractionModifiers coulomb_modifier = InteractionModifiers::Default;
+ //! Coulomb switch range start (nm)
+ real rcoulomb_switch = 0;
+ //! Coulomb cutoff (nm)
+ real rcoulomb = 0;
+ //! Relative dielectric constant
+ real epsilon_r = 0;
+ //! Relative dielectric constant of the RF
+ real epsilon_rf = 0;
+ //! Always false (no longer supported)
+ bool implicit_solvent = false;
+ //! Type of Van der Waals treatment
+ VanDerWaalsType vdwtype = VanDerWaalsType::Default;
+ //! Modify the Van der Waals interaction
+ InteractionModifiers vdw_modifier = InteractionModifiers::Default;
+ //! Van der Waals switch range start (nm)
+ real rvdw_switch = 0;
+ //! Van der Waals cutoff (nm)
+ real rvdw = 0;
+ //! Perform Long range dispersion corrections
+ DispersionCorrectionType eDispCorr = DispersionCorrectionType::Default;
+ //! Extension of the table beyond the cut-off, as well as the table length for 1-4 interac.
+ real tabext = 0;
+ //! Tolerance for shake
+ real shake_tol = 0;
+ //! Free energy calculations
+ FreeEnergyPerturbationType efep = FreeEnergyPerturbationType::Default;
+ //! Data for the FEP state
+ std::unique_ptr<t_lambda> fepvals;
+ //! Whether to do simulated tempering
+ bool bSimTemp = false;
+ //! Variables for simulated tempering
+ std::unique_ptr<t_simtemp> simtempvals;
+ //! Whether expanded ensembles are used
+ bool bExpanded = false;
+ //! Expanded ensemble parameters
+ std::unique_ptr<t_expanded> expandedvals;
+ //! Type of distance restraining
+ DistanceRestraintRefinement eDisre = DistanceRestraintRefinement::Default;
+ //! Force constant for time averaged distance restraints
+ real dr_fc = 0;
+ //! Type of weighting of pairs in one restraints
+ DistanceRestraintWeighting eDisreWeighting = DistanceRestraintWeighting::Default;
+ //! Use combination of time averaged and instantaneous violations
+ bool bDisreMixed = false;
+ //! Frequency of writing pair distances to enx
+ int nstdisreout = 0;
+ //! Time constant for memory function in disres
+ real dr_tau = 0;
+ //! Force constant for orientational restraints
+ real orires_fc = 0;
+ //! Time constant for memory function in orires
+ real orires_tau = 0;
+ //! Frequency of writing tr(SD) to energy output
+ int nstorireout = 0;
+ //! The stepsize for updating
+ real em_stepsize = 0;
+ //! The tolerance
+ real em_tol = 0;
+ //! Number of iterations for convergence of steepest descent in relax_shells
+ int niter = 0;
+ //! Stepsize for directional minimization in relax_shells
+ real fc_stepsize = 0;
+ //! Number of steps after which a steepest descents step is done while doing cg
+ int nstcgsteep = 0;
+ //! Number of corrections to the Hessian to keep
+ int nbfgscorr = 0;
+ //! Type of constraint algorithm
+ ConstraintAlgorithm eConstrAlg = ConstraintAlgorithm::Default;
+ //! Order of the LINCS Projection Algorithm
+ int nProjOrder = 0;
+ //! Warn if any bond rotates more than this many degrees
+ real LincsWarnAngle = 0;
+ //! Number of iterations in the final LINCS step
+ int nLincsIter = 0;
+ //! Use successive overrelaxation for shake
+ bool bShakeSOR = false;
+ //! Friction coefficient for BD (amu/ps)
+ real bd_fric = 0;
+ //! Random seed for SD and BD
+ int64_t ld_seed = 0;
+ //! The number of walls
+ int nwall = 0;
+ //! The type of walls
+ WallType wall_type = WallType::Default;
+ //! The potentail is linear for r<=wall_r_linpot
+ real wall_r_linpot = 0;
+ //! The atom type for walls
+ int wall_atomtype[2] = { 0, 0 };
+ //! Number density for walls
+ real wall_density[2] = { 0, 0 };
+ //! Scaling factor for the box for Ewald
+ real wall_ewald_zfac = 0;
/* COM pulling data */
- gmx_bool bPull; /* Do we do COM pulling? */
- struct pull_params_t *pull; /* The data for center of mass pulling */
- struct pull_t *pull_work; /* The COM pull force calculation data structure; TODO this pointer should live somewhere else */
+ //! Do we do COM pulling?
+ bool bPull = false;
+ //! The data for center of mass pulling
+ std::unique_ptr<pull_params_t> pull;
+
+ /* AWH bias data */
+ //! Whether to use AWH biasing for PMF calculations
+ bool bDoAwh = false;
+ //! AWH biasing parameters
+ std::unique_ptr<gmx::AwhParams> awhParams;
/* Enforced rotation data */
- gmx_bool bRot; /* Calculate enforced rotation potential(s)? */
- t_rot *rot; /* The data for enforced rotation potentials */
-
- int eSwapCoords; /* Do ion/water position exchanges (CompEL)? */
- t_swapcoords *swap;
-
- gmx_bool bIMD; /* Allow interactive MD sessions for this .tpr? */
- t_IMD *imd; /* Interactive molecular dynamics */
-
- real cos_accel; /* Acceleration for viscosity calculation */
- tensor deform; /* Triclinic deformation velocities (nm/ps) */
- int userint1; /* User determined parameters */
- int userint2;
- int userint3;
- int userint4;
- real userreal1;
- real userreal2;
- real userreal3;
- real userreal4;
- t_grpopts opts; /* Group options */
- gmx_bool bQMMM; /* QM/MM calculation */
- int QMconstraints; /* constraints on QM bonds */
- int QMMMscheme; /* Scheme: ONIOM or normal */
- real scalefactor; /* factor for scaling the MM charges in QM calc.*/
+ //! Whether to calculate enforced rotation potential(s)
+ bool bRot = false;
+ //! The data for enforced rotation potentials
+ t_rot* rot = nullptr;
+
+ //! Whether to do ion/water position exchanges (CompEL)
+ SwapType eSwapCoords = SwapType::Default;
+ //! Swap data structure.
+ t_swapcoords* swap = nullptr;
+
+ //! Whether the tpr makes an interactive MD session possible.
+ bool bIMD = false;
+ //! Interactive molecular dynamics
+ t_IMD* imd = nullptr;
+
+ //! Acceleration for viscosity calculation
+ real cos_accel = 0;
+ //! Triclinic deformation velocities (nm/ps)
+ tensor deform = { { 0 } };
+ /*! \brief User determined parameters */
+ /**@{*/
+ int userint1 = 0;
+ int userint2 = 0;
+ int userint3 = 0;
+ int userint4 = 0;
+ real userreal1 = 0;
+ real userreal2 = 0;
+ real userreal3 = 0;
+ real userreal4 = 0;
+ /**@}*/
+ //! Group options
+ t_grpopts opts;
+ //! QM/MM calculation
+ bool bQMMM = false;
/* Fields for removed features go here (better caching) */
- gmx_bool bAdress; // Whether AdResS is enabled - always false if a valid .tpr was read
- gmx_bool useTwinRange; // Whether twin-range scheme is active - always false if a valid .tpr was read
-
- gmx::KeyValueTreeObject *params;
+ //! Whether AdResS is enabled - always false if a valid .tpr was read
+ bool bAdress = false;
+ //! Whether twin-range scheme is active - always false if a valid .tpr was read
+ bool useTwinRange = false;
+ //! Whether we have constant acceleration
+ bool useConstantAcceleration = false;
+
+ //! KVT object that contains input parameters converted to the new style.
+ gmx::KeyValueTreeObject* params = nullptr;
+
+ //! KVT for storing simulation parameters that are not part of the mdp file.
+ std::unique_ptr<gmx::KeyValueTreeObject> internalParameters;
};
-int ir_optimal_nstcalcenergy(const t_inputrec *ir);
+int ir_optimal_nstcalcenergy(const t_inputrec* ir);
-int tcouple_min_integration_steps(int etc);
+int tcouple_min_integration_steps(TemperatureCoupling etc);
-int ir_optimal_nsttcouple(const t_inputrec *ir);
+int ir_optimal_nsttcouple(const t_inputrec* ir);
-int pcouple_min_integration_steps(int epc);
+int pcouple_min_integration_steps(PressureCoupling epc);
-int ir_optimal_nstpcouple(const t_inputrec *ir);
+int ir_optimal_nstpcouple(const t_inputrec* ir);
/* Returns if the Coulomb force or potential is switched to zero */
-gmx_bool ir_coulomb_switched(const t_inputrec *ir);
+bool ir_coulomb_switched(const t_inputrec* ir);
/* Returns if the Coulomb interactions are zero beyond the rcoulomb.
* Note: always returns TRUE for the Verlet cut-off scheme.
*/
-gmx_bool ir_coulomb_is_zero_at_cutoff(const t_inputrec *ir);
+bool ir_coulomb_is_zero_at_cutoff(const t_inputrec* ir);
/* As ir_coulomb_is_zero_at_cutoff, but also returns TRUE for user tabulated
* interactions, since these might be zero beyond rcoulomb.
*/
-gmx_bool ir_coulomb_might_be_zero_at_cutoff(const t_inputrec *ir);
+bool ir_coulomb_might_be_zero_at_cutoff(const t_inputrec* ir);
/* Returns if the Van der Waals force or potential is switched to zero */
-gmx_bool ir_vdw_switched(const t_inputrec *ir);
+bool ir_vdw_switched(const t_inputrec* ir);
/* Returns if the Van der Waals interactions are zero beyond the rvdw.
* Note: always returns TRUE for the Verlet cut-off scheme.
*/
-gmx_bool ir_vdw_is_zero_at_cutoff(const t_inputrec *ir);
+bool ir_vdw_is_zero_at_cutoff(const t_inputrec* ir);
/* As ir_vdw_is_zero_at_cutoff, but also returns TRUE for user tabulated
* interactions, since these might be zero beyond rvdw.
*/
-gmx_bool ir_vdw_might_be_zero_at_cutoff(const t_inputrec *ir);
+bool ir_vdw_might_be_zero_at_cutoff(const t_inputrec* ir);
/*! \brief Free memory from input record.
*
*
* \param[in] ir The data structure
*/
-void done_inputrec(t_inputrec *ir);
+void done_inputrec(t_inputrec* ir);
+
+void pr_inputrec(FILE* fp, int indent, const char* title, const t_inputrec* ir, bool bMDPformat);
-void pr_inputrec(FILE *fp, int indent, const char *title, const t_inputrec *ir,
- gmx_bool bMDPformat);
+void cmp_inputrec(FILE* fp, const t_inputrec* ir1, const t_inputrec* ir2, real ftol, real abstol);
-void cmp_inputrec(FILE *fp, const t_inputrec *ir1, const t_inputrec *ir2, real ftol, real abstol);
+void comp_pull_AB(FILE* fp, const pull_params_t& pull, real ftol, real abstol);
-void comp_pull_AB(FILE *fp, pull_params_t *pull, real ftol, real abstol);
+bool inputrecDeform(const t_inputrec* ir);
-gmx_bool inputrecDeform(const t_inputrec *ir);
+bool inputrecDynamicBox(const t_inputrec* ir);
-gmx_bool inputrecDynamicBox(const t_inputrec *ir);
+bool inputrecPreserveShape(const t_inputrec* ir);
-gmx_bool inputrecPreserveShape(const t_inputrec *ir);
+bool inputrecNeedMutot(const t_inputrec* ir);
-gmx_bool inputrecNeedMutot(const t_inputrec *ir);
+bool inputrecTwinRange(const t_inputrec* ir);
-gmx_bool inputrecTwinRange(const t_inputrec *ir);
+bool inputrecExclForces(const t_inputrec* ir);
-gmx_bool inputrecExclForces(const t_inputrec *ir);
+bool inputrecNptTrotter(const t_inputrec* ir);
-gmx_bool inputrecNptTrotter(const t_inputrec *ir);
+bool inputrecNvtTrotter(const t_inputrec* ir);
-gmx_bool inputrecNvtTrotter(const t_inputrec *ir);
+bool inputrecNphTrotter(const t_inputrec* ir);
-gmx_bool inputrecNphTrotter(const t_inputrec *ir);
+/*! \brief Return true if the simulation is 2D periodic with two walls. */
+bool inputrecPbcXY2Walls(const t_inputrec* ir);
-/* Returns true for MD integator with T and/or P-coupling that supports
- * calculating the conserved energy quantity.
+//! \brief Return true if the simulation has frozen atoms (non-trivial freeze groups).
+bool inputrecFrozenAtoms(const t_inputrec* ir);
+
+/*! \brief Returns true for MD integator with T and/or P-coupling that supports
+ * calculating a conserved energy quantity.
+ *
+ * Note that dynamical integrators without T and P coupling (ie NVE)
+ * return false, i.e. the return value refers to whether there
+ * is a conserved quantity different than the total energy.
*/
-bool integratorHasConservedEnergyQuantity(const t_inputrec *ir);
+bool integratorHasConservedEnergyQuantity(const t_inputrec* ir);
+
+/*! \brief Returns true when temperature is coupled or constant. */
+bool integratorHasReferenceTemperature(const t_inputrec* ir);
/*! \brief Return the number of bounded dimensions
*
* \return the number of dimensions in which
* the coordinates of the particles are bounded, starting at X.
*/
-int inputrec2nboundeddim(const t_inputrec *ir);
+int inputrec2nboundeddim(const t_inputrec* ir);
/*! \brief Returns the number of degrees of freedom in center of mass motion
*
- * \param[in] ir the inputrec structure
+ * \param[in] ir The inputrec structure
* \return the number of degrees of freedom of the center of mass
*/
-int ndof_com(const t_inputrec *ir);
+int ndof_com(const t_inputrec* ir);
+
+/*! \brief Returns the maximum reference temperature over all coupled groups
+ *
+ * Returns 0 for energy minimization and normal mode computation.
+ * Returns -1 for MD without temperature coupling.
+ *
+ * \param[in] ir The inputrec structure
+ */
+real maxReferenceTemperature(const t_inputrec& ir);
+
+/*! \brief Returns whether there is an Ewald surface contribution
+ */
+bool haveEwaldSurfaceContribution(const t_inputrec& ir);
+
+/*! \brief Check if calculation of the specific FEP type was requested.
+ *
+ * \param[in] ir Input record.
+ * \param[in] fepType Free-energy perturbation type to check for.
+ *
+ * \returns If the \p fepType is perturbed in this run.
+ */
+bool haveFreeEnergyType(const t_inputrec& ir, int fepType);
#endif /* GMX_MDTYPES_INPUTREC_H */
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff