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, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
38 * Declares inputrec data structure and utilities.
41 * \ingroup module_mdtypes
43 #ifndef GMX_MDTYPES_INPUTREC_H
44 #define GMX_MDTYPES_INPUTREC_H
48 #include "gromacs/math/vectypes.h"
49 #include "gromacs/mdtypes/md_enums.h"
50 #include "gromacs/utility/basedefinitions.h"
51 #include "gromacs/utility/real.h"
56 //! Coeffients (V / nm)
63 real E0
; /* Field strength (V/nm) */
64 real omega
; /* Frequency (1/ps) */
65 real t0
; /* Centre of the Gaussian pulse (ps) */
66 real sigma
; /* Width of the Gaussian pulse (FWHM) (ps) */
69 #define EGP_EXCL (1<<0)
70 #define EGP_TABLE (1<<1)
72 typedef struct t_grpopts
{
73 int ngtc
; /* # T-Coupl groups */
74 int nhchainlength
; /* # of nose-hoover chains per group */
75 int ngacc
; /* # Accelerate groups */
76 int ngfrz
; /* # Freeze groups */
77 int ngener
; /* # Ener groups */
78 real
*nrdf
; /* Nr of degrees of freedom in a group */
79 real
*ref_t
; /* Coupling temperature per group */
80 int *annealing
; /* No/simple/periodic SA for each group */
81 int *anneal_npoints
; /* Number of annealing time points per grp */
82 real
**anneal_time
; /* For ea. group: Time points */
83 real
**anneal_temp
; /* For ea. grp: Temperature at these times */
84 /* Final temp after all intervals is ref_t */
85 real
*tau_t
; /* Tau coupling time */
86 rvec
*acc
; /* Acceleration per group */
87 ivec
*nFreeze
; /* Freeze the group in each direction ? */
88 int *egp_flags
; /* Exclusions/tables of energy group pairs */
91 int ngQM
; /* nr of QM groups */
92 int *QMmethod
; /* Level of theory in the QM calculation */
93 int *QMbasis
; /* Basisset in the QM calculation */
94 int *QMcharge
; /* Total charge in the QM region */
95 int *QMmult
; /* Spin multiplicicty in the QM region */
96 gmx_bool
*bSH
; /* surface hopping (diabatic hop only) */
97 int *CASorbitals
; /* number of orbiatls in the active space */
98 int *CASelectrons
; /* number of electrons in the active space */
99 real
*SAon
; /* at which gap (A.U.) the SA is switched on */
101 int *SAsteps
; /* in how many steps SA goes from 1-1 to 0.5-0.5*/
106 typedef struct t_simtemp
{
107 int eSimTempScale
; /* simulated temperature scaling; linear or exponential */
108 real simtemp_low
; /* the low temperature for simulated tempering */
109 real simtemp_high
; /* the high temperature for simulated tempering */
110 real
*temperatures
; /* the range of temperatures used for simulated tempering */
113 typedef struct t_lambda
{
114 int nstdhdl
; /* The frequency for calculating dhdl */
115 double init_lambda
; /* fractional value of lambda (usually will use
116 init_fep_state, this will only be for slow growth,
117 and for legacy free energy code. Only has a
118 valid value if positive) */
119 int init_fep_state
; /* the initial number of the state */
120 double delta_lambda
; /* change of lambda per time step (fraction of (0.1) */
121 int edHdLPrintEnergy
; /* print no, total or potential energies in dhdl */
122 int n_lambda
; /* The number of foreign lambda points */
123 double **all_lambda
; /* The array of all lambda values */
124 int lambda_neighbors
; /* The number of neighboring lambda states to
125 calculate the energy for in up and down directions
127 int lambda_start_n
; /* The first lambda to calculate energies for */
128 int lambda_stop_n
; /* The last lambda +1 to calculate energies for */
129 real sc_alpha
; /* free energy soft-core parameter */
130 int sc_power
; /* lambda power for soft-core interactions */
131 real sc_r_power
; /* r power for soft-core interactions */
132 real sc_sigma
; /* free energy soft-core sigma when c6 or c12=0 */
133 real sc_sigma_min
; /* free energy soft-core sigma for ????? */
134 gmx_bool bScCoul
; /* use softcore for the coulomb portion as well (default FALSE) */
135 gmx_bool separate_dvdl
[efptNR
]; /* whether to print the dvdl term associated with
136 this term; if it is not specified as separate,
137 it is lumped with the FEP term */
138 int separate_dhdl_file
; /* whether to write a separate dhdl.xvg file
139 note: NOT a gmx_bool, but an enum */
140 int dhdl_derivatives
; /* whether to calculate+write dhdl derivatives
141 note: NOT a gmx_bool, but an enum */
142 int dh_hist_size
; /* The maximum table size for the dH histogram */
143 double dh_hist_spacing
; /* The spacing for the dH histogram */
146 typedef struct t_expanded
{
147 int nstexpanded
; /* The frequency of expanded ensemble state changes */
148 int elamstats
; /* which type of move updating do we use for lambda monte carlo (or no for none) */
149 int elmcmove
; /* what move set will be we using for state space moves */
150 int elmceq
; /* the method we use to decide of we have equilibrated the weights */
151 int equil_n_at_lam
; /* the minumum number of samples at each lambda for deciding whether we have reached a minimum */
152 real equil_wl_delta
; /* WL delta at which we stop equilibrating weights */
153 real equil_ratio
; /* use the ratio of weights (ratio of minimum to maximum) to decide when to stop equilibrating */
154 int equil_steps
; /* after equil_steps steps we stop equilibrating the weights */
155 int equil_samples
; /* after equil_samples total samples (steps/nstfep), we stop equilibrating the weights */
156 int lmc_seed
; /* random number seed for lambda mc switches */
157 gmx_bool minvar
; /* whether to use minumum variance weighting */
158 int minvarmin
; /* the number of samples needed before kicking into minvar routine */
159 real minvar_const
; /* the offset for the variance in MinVar */
160 int c_range
; /* range of cvalues used for BAR */
161 gmx_bool bSymmetrizedTMatrix
; /* whether to print symmetrized matrices */
162 int nstTij
; /* How frequently to print the transition matrices */
163 int lmc_repeats
; /* number of repetitions in the MC lambda jumps */ /*MRS -- VERIFY THIS */
164 int lmc_forced_nstart
; /* minimum number of samples for each state before free sampling */ /* MRS -- VERIFY THIS! */
165 int gibbsdeltalam
; /* distance in lambda space for the gibbs interval */
166 real wl_scale
; /* scaling factor for wang-landau */
167 real wl_ratio
; /* ratio between largest and smallest number for freezing the weights */
168 real init_wl_delta
; /* starting delta for wang-landau */
169 gmx_bool bWLoneovert
; /* use one over t convergece for wang-landau when the delta get sufficiently small */
170 gmx_bool bInit_weights
; /* did we initialize the weights? TODO: REMOVE FOR 5.0, no longer needed with new logic */
171 real mc_temp
; /* To override the main temperature, or define it if it's not defined */
172 real
*init_lambda_weights
; /* user-specified initial weights to start with */
176 /* Abstract types for enforced rotation only defined in pull_rotation.c */
177 typedef struct gmx_enfrot
*gmx_enfrot_t
;
178 typedef struct gmx_enfrotgrp
*gmx_enfrotgrp_t
;
181 int eType
; /* Rotation type for this group */
182 int bMassW
; /* Use mass-weighed positions? */
183 int nat
; /* Number of atoms in the group */
184 int *ind
; /* The global atoms numbers */
185 rvec
*x_ref
; /* The reference positions */
186 rvec vec
; /* The normalized rotation vector */
187 real rate
; /* Rate of rotation (degree/ps) */
188 real k
; /* Force constant (kJ/(mol nm^2) */
189 rvec pivot
; /* Pivot point of rotation axis (nm) */
190 int eFittype
; /* Type of fit to determine actual group angle */
191 int PotAngle_nstep
; /* Number of angles around the reference angle
192 for which the rotation potential is also
193 evaluated (for fit type 'potential' only) */
194 real PotAngle_step
; /* Distance between two angles in degrees (for
195 fit type 'potential' only) */
196 real slab_dist
; /* Slab distance (nm) */
197 real min_gaussian
; /* Minimum value the gaussian must have so that
198 the force is actually evaluated */
199 real eps
; /* Additive constant for radial motion2 and
200 flexible2 potentials (nm^2) */
201 gmx_enfrotgrp_t enfrotgrp
; /* Stores non-inputrec rotation data per group */
204 typedef struct t_rot
{
205 int ngrp
; /* Number of rotation groups */
206 int nstrout
; /* Output frequency for main rotation outfile */
207 int nstsout
; /* Output frequency for per-slab data */
208 t_rotgrp
*grp
; /* Groups to rotate */
209 gmx_enfrot_t enfrot
; /* Stores non-inputrec enforced rotation data */
212 /* Abstract type for IMD only defined in IMD.c */
215 typedef struct t_IMD
{
216 int nat
; /* Number of interactive atoms */
217 int *ind
; /* The global indices of the interactive atoms */
218 struct t_gmx_IMD
*setup
; /* Stores non-inputrec IMD data */
221 /* Abstract types for position swapping only defined in swapcoords.cpp */
222 typedef struct t_swap
*gmx_swapcoords_t
;
224 typedef struct t_swapGroup
{
225 char *molname
; /* Name of the swap group, e.g. NA, CL, SOL */
226 int nat
; /* Number of atoms in this group */
227 int *ind
; /* The global ion group atoms numbers */
228 int nmolReq
[eCompNR
]; /* Requested number of molecules of this type
232 typedef struct t_swapcoords
{
233 int nstswap
; /* Every how many steps a swap is attempted? */
234 gmx_bool massw_split
[2]; /* Use mass-weighted positions in split group? */
235 real cyl0r
, cyl1r
; /* Split cylinders defined by radius, upper and */
236 real cyl0u
, cyl1u
; /* ... lower extension. The split cylinders de- */
237 real cyl0l
, cyl1l
; /* ... fine the channels and are each anchored */
238 /* ... in the center of the split group */
239 int nAverage
; /* Coupling constant (nr of swap attempt steps) */
240 real threshold
; /* Ion counts may deviate from the requested
241 values by +-threshold before a swap is done */
242 real bulkOffset
[eCompNR
]; /* Offset of the swap layer (='bulk') w.r.t.
243 the compartment-defining layers */
244 int ngrp
; /* Number of groups to be controlled */
245 t_swapGroup
*grp
; /* All swap groups, including split and solvent */
246 gmx_swapcoords_t si_priv
; /* swap private data accessible in
250 typedef struct t_inputrec
{
251 int eI
; /* Integration method */
252 gmx_int64_t nsteps
; /* number of steps to be taken */
253 int simulation_part
; /* Used in checkpointing to separate chunks */
254 gmx_int64_t init_step
; /* start at a stepcount >0 (used w. convert-tpr) */
255 int nstcalcenergy
; /* frequency of energy calc. and T/P coupl. upd. */
256 int cutoff_scheme
; /* group or verlet cutoffs */
257 int ns_type
; /* which ns method should we use? */
258 int nstlist
; /* number of steps before pairlist is generated */
259 int ndelta
; /* number of cells per rlong */
260 int nstcomm
; /* number of steps after which center of mass */
261 /* motion is removed */
262 int comm_mode
; /* Center of mass motion removal algorithm */
263 int nstlog
; /* number of steps after which print to logfile */
264 int nstxout
; /* number of steps after which X is output */
265 int nstvout
; /* id. for V */
266 int nstfout
; /* id. for F */
267 int nstenergy
; /* number of steps after which energies printed */
268 int nstxout_compressed
; /* id. for compressed trj (.xtc,.tng) */
269 double init_t
; /* initial time (ps) */
270 double delta_t
; /* time step (ps) */
271 real x_compression_precision
; /* precision of x in compressed trajectory file */
272 real fourier_spacing
; /* requested fourier_spacing, when nk? not set */
273 int nkx
, nky
, nkz
; /* number of k vectors in each spatial dimension*/
274 /* for fourier methods for long range electrost.*/
275 int pme_order
; /* interpolation order for PME */
276 real ewald_rtol
; /* Real space tolerance for Ewald, determines */
277 /* the real/reciprocal space relative weight */
278 real ewald_rtol_lj
; /* Real space tolerance for LJ-Ewald */
279 int ewald_geometry
; /* normal/3d ewald, or pseudo-2d LR corrections */
280 real epsilon_surface
; /* Epsilon for PME dipole correction */
281 int ljpme_combination_rule
; /* Type of combination rule in LJ-PME */
282 int ePBC
; /* Type of periodic boundary conditions */
283 int bPeriodicMols
; /* Periodic molecules */
284 gmx_bool bContinuation
; /* Continuation run: starting state is correct */
285 int etc
; /* temperature coupling */
286 int nsttcouple
; /* interval in steps for temperature coupling */
287 gmx_bool bPrintNHChains
; /* whether to print nose-hoover chains */
288 int epc
; /* pressure coupling */
289 int epct
; /* pressure coupling type */
290 int nstpcouple
; /* interval in steps for pressure coupling */
291 real tau_p
; /* pressure coupling time (ps) */
292 tensor ref_p
; /* reference pressure (kJ/(mol nm^3)) */
293 tensor compress
; /* compressability ((mol nm^3)/kJ) */
294 int refcoord_scaling
; /* How to scale absolute reference coordinates */
295 rvec posres_com
; /* The COM of the posres atoms */
296 rvec posres_comB
; /* The B-state COM of the posres atoms */
297 int andersen_seed
; /* Random seed for Andersen thermostat (obsolete) */
298 real verletbuf_tol
; /* Per atom pair energy drift tolerance (kJ/mol/ps/atom) for list buffer */
299 real rlist
; /* short range pairlist cut-off (nm) */
300 real rtpi
; /* Radius for test particle insertion */
301 int coulombtype
; /* Type of electrostatics treatment */
302 int coulomb_modifier
; /* Modify the Coulomb interaction */
303 real rcoulomb_switch
; /* Coulomb switch range start (nm) */
304 real rcoulomb
; /* Coulomb cutoff (nm) */
305 real epsilon_r
; /* relative dielectric constant */
306 real epsilon_rf
; /* relative dielectric constant of the RF */
307 int implicit_solvent
; /* No (=explicit water), or GBSA solvent models */
308 int gb_algorithm
; /* Algorithm to use for calculation Born radii */
309 int nstgbradii
; /* Frequency of updating Generalized Born radii */
310 real rgbradii
; /* Cutoff for GB radii calculation */
311 real gb_saltconc
; /* Salt concentration (M) for GBSA models */
312 real gb_epsilon_solvent
; /* dielectric coeff. of implicit solvent */
313 real gb_obc_alpha
; /* 1st scaling factor for Bashford-Case GB */
314 real gb_obc_beta
; /* 2nd scaling factor for Bashford-Case GB */
315 real gb_obc_gamma
; /* 3rd scaling factor for Bashford-Case GB */
316 real gb_dielectric_offset
; /* Dielectric offset for Still/HCT/OBC */
317 int sa_algorithm
; /* Algorithm for SA part of GBSA */
318 real sa_surface_tension
; /* Energy factor for SA part of GBSA */
319 int vdwtype
; /* Type of Van der Waals treatment */
320 int vdw_modifier
; /* Modify the VdW interaction */
321 real rvdw_switch
; /* Van der Waals switch range start (nm) */
322 real rvdw
; /* Van der Waals cutoff (nm) */
323 int eDispCorr
; /* Perform Long range dispersion corrections */
324 real tabext
; /* Extension of the table beyond the cut-off, *
325 * as well as the table length for 1-4 interac. */
326 real shake_tol
; /* tolerance for shake */
327 int efep
; /* free energy calculations */
328 t_lambda
*fepvals
; /* Data for the FEP state */
329 gmx_bool bSimTemp
; /* Whether to do simulated tempering */
330 t_simtemp
*simtempvals
; /* Variables for simulated tempering */
331 gmx_bool bExpanded
; /* Whether expanded ensembles are used */
332 t_expanded
*expandedvals
; /* Expanded ensemble parameters */
333 int eDisre
; /* Type of distance restraining */
334 real dr_fc
; /* force constant for ta_disre */
335 int eDisreWeighting
; /* type of weighting of pairs in one restraints */
336 gmx_bool bDisreMixed
; /* Use comb of time averaged and instan. viol's */
337 int nstdisreout
; /* frequency of writing pair distances to enx */
338 real dr_tau
; /* time constant for memory function in disres */
339 real orires_fc
; /* force constant for orientational restraints */
340 real orires_tau
; /* time constant for memory function in orires */
341 int nstorireout
; /* frequency of writing tr(SD) to enx */
342 real em_stepsize
; /* The stepsize for updating */
343 real em_tol
; /* The tolerance */
344 int niter
; /* Number of iterations for convergence of */
345 /* steepest descent in relax_shells */
346 real fc_stepsize
; /* Stepsize for directional minimization */
347 /* in relax_shells */
348 int nstcgsteep
; /* number of steps after which a steepest */
349 /* descents step is done while doing cg */
350 int nbfgscorr
; /* Number of corrections to the hessian to keep */
351 int eConstrAlg
; /* Type of constraint algorithm */
352 int nProjOrder
; /* Order of the LINCS Projection Algorithm */
353 real LincsWarnAngle
; /* If bond rotates more than %g degrees, warn */
354 int nLincsIter
; /* Number of iterations in the final Lincs step */
355 gmx_bool bShakeSOR
; /* Use successive overrelaxation for shake */
356 real bd_fric
; /* Friction coefficient for BD (amu/ps) */
357 gmx_int64_t ld_seed
; /* Random seed for SD and BD */
358 int nwall
; /* The number of walls */
359 int wall_type
; /* The type of walls */
360 real wall_r_linpot
; /* The potentail is linear for r<=wall_r_linpot */
361 int wall_atomtype
[2]; /* The atom type for walls */
362 real wall_density
[2]; /* Number density for walls */
363 real wall_ewald_zfac
; /* Scaling factor for the box for Ewald */
365 /* COM pulling data */
366 gmx_bool bPull
; /* Do we do COM pulling? */
367 struct pull_params_t
*pull
; /* The data for center of mass pulling */
368 struct pull_t
*pull_work
; /* The COM pull force calculation data structure; TODO this pointer should live somewhere else */
370 /* Enforced rotation data */
371 gmx_bool bRot
; /* Calculate enforced rotation potential(s)? */
372 t_rot
*rot
; /* The data for enforced rotation potentials */
374 int eSwapCoords
; /* Do ion/water position exchanges (CompEL)? */
377 gmx_bool bIMD
; /* Allow interactive MD sessions for this .tpr? */
378 t_IMD
*imd
; /* Interactive molecular dynamics */
380 real cos_accel
; /* Acceleration for viscosity calculation */
381 tensor deform
; /* Triclinic deformation velocities (nm/ps) */
382 int userint1
; /* User determined parameters */
390 t_grpopts opts
; /* Group options */
391 t_cosines ex
[DIM
]; /* Electric field stuff (spatial part) */
392 t_cosines et
[DIM
]; /* Electric field stuff (time part) */
393 gmx_bool bQMMM
; /* QM/MM calculation */
394 int QMconstraints
; /* constraints on QM bonds */
395 int QMMMscheme
; /* Scheme: ONIOM or normal */
396 real scalefactor
; /* factor for scaling the MM charges in QM calc.*/
398 /* Fields for removed features go here (better caching) */
399 gmx_bool bAdress
; // Whether AdResS is enabled - always false if a valid .tpr was read
400 gmx_bool useTwinRange
; // Whether twin-range scheme is active - always false if a valid .tpr was read
403 int ir_optimal_nstcalcenergy(const t_inputrec
*ir
);
405 int tcouple_min_integration_steps(int etc
);
407 int ir_optimal_nsttcouple(const t_inputrec
*ir
);
409 int pcouple_min_integration_steps(int epc
);
411 int ir_optimal_nstpcouple(const t_inputrec
*ir
);
413 /* Returns if the Coulomb force or potential is switched to zero */
414 gmx_bool
ir_coulomb_switched(const t_inputrec
*ir
);
416 /* Returns if the Coulomb interactions are zero beyond the rcoulomb.
417 * Note: always returns TRUE for the Verlet cut-off scheme.
419 gmx_bool
ir_coulomb_is_zero_at_cutoff(const t_inputrec
*ir
);
421 /* As ir_coulomb_is_zero_at_cutoff, but also returns TRUE for user tabulated
422 * interactions, since these might be zero beyond rcoulomb.
424 gmx_bool
ir_coulomb_might_be_zero_at_cutoff(const t_inputrec
*ir
);
426 /* Returns if the Van der Waals force or potential is switched to zero */
427 gmx_bool
ir_vdw_switched(const t_inputrec
*ir
);
429 /* Returns if the Van der Waals interactions are zero beyond the rvdw.
430 * Note: always returns TRUE for the Verlet cut-off scheme.
432 gmx_bool
ir_vdw_is_zero_at_cutoff(const t_inputrec
*ir
);
434 /* As ir_vdw_is_zero_at_cutoff, but also returns TRUE for user tabulated
435 * interactions, since these might be zero beyond rvdw.
437 gmx_bool
ir_vdw_might_be_zero_at_cutoff(const t_inputrec
*ir
);
439 /*! \brief Initiate input record structure
441 * Initialiazes all the arrays and pointers to NULL.
443 * \param[in] ir Inputrec must be pre-allocated
445 void init_inputrec(t_inputrec
*ir
);
447 /*! \brief Free memory from input record.
449 * All arrays and pointers will be freed.
451 * \param[in] ir The data structure
453 void done_inputrec(t_inputrec
*ir
);
455 void pr_inputrec(FILE *fp
, int indent
, const char *title
, const t_inputrec
*ir
,
456 gmx_bool bMDPformat
);
458 gmx_bool
inputrecDeform(const t_inputrec
*ir
);
460 gmx_bool
inputrecDynamicBox(const t_inputrec
*ir
);
462 gmx_bool
inputrecPreserveShape(const t_inputrec
*ir
);
464 gmx_bool
inputrecNeedMutot(const t_inputrec
*ir
);
466 gmx_bool
inputrecTwinRange(const t_inputrec
*ir
);
468 gmx_bool
inputrecElecField(const t_inputrec
*ir
);
470 gmx_bool
inputrecExclForces(const t_inputrec
*ir
);
472 gmx_bool
inputrecNptTrotter(const t_inputrec
*ir
);
474 gmx_bool
inputrecNvtTrotter(const t_inputrec
*ir
);
476 gmx_bool
inputrecNphTrotter(const t_inputrec
*ir
);
478 #endif /* GMX_MDTYPES_INPUTREC_H */
