Update Awh initialization and lifetime management
[gromacs.git] / src / gromacs / mdtypes / inputrec.h
blob3eff130d5531577d69ed86e5fa296a23a7c0eb07
1 /*
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,2018, 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.
37 #ifndef GMX_MDTYPES_INPUTREC_H
38 #define GMX_MDTYPES_INPUTREC_H
40 #include <cstdio>
42 #include <memory>
44 #include "gromacs/math/vectypes.h"
45 #include "gromacs/mdtypes/md_enums.h"
46 #include "gromacs/utility/basedefinitions.h"
47 #include "gromacs/utility/real.h"
49 #define EGP_EXCL (1<<0)
50 #define EGP_TABLE (1<<1)
52 struct pull_params_t;
54 namespace gmx
56 class Awh;
57 struct AwhParams;
58 class KeyValueTreeObject;
61 typedef struct t_grpopts {
62 int ngtc; /* # T-Coupl groups */
63 int nhchainlength; /* # of nose-hoover chains per group */
64 int ngacc; /* # Accelerate groups */
65 int ngfrz; /* # Freeze groups */
66 int ngener; /* # Ener groups */
67 real *nrdf; /* Nr of degrees of freedom in a group */
68 real *ref_t; /* Coupling temperature per group */
69 int *annealing; /* No/simple/periodic SA for each group */
70 int *anneal_npoints; /* Number of annealing time points per grp */
71 real **anneal_time; /* For ea. group: Time points */
72 real **anneal_temp; /* For ea. grp: Temperature at these times */
73 /* Final temp after all intervals is ref_t */
74 real *tau_t; /* Tau coupling time */
75 rvec *acc; /* Acceleration per group */
76 ivec *nFreeze; /* Freeze the group in each direction ? */
77 int *egp_flags; /* Exclusions/tables of energy group pairs */
79 /* QMMM stuff */
80 int ngQM; /* nr of QM groups */
81 int *QMmethod; /* Level of theory in the QM calculation */
82 int *QMbasis; /* Basisset in the QM calculation */
83 int *QMcharge; /* Total charge in the QM region */
84 int *QMmult; /* Spin multiplicicty in the QM region */
85 gmx_bool *bSH; /* surface hopping (diabatic hop only) */
86 int *CASorbitals; /* number of orbiatls in the active space */
87 int *CASelectrons; /* number of electrons in the active space */
88 real *SAon; /* at which gap (A.U.) the SA is switched on */
89 real *SAoff;
90 int *SAsteps; /* in how many steps SA goes from 1-1 to 0.5-0.5*/
91 } t_grpopts;
93 typedef struct t_simtemp {
94 int eSimTempScale; /* simulated temperature scaling; linear or exponential */
95 real simtemp_low; /* the low temperature for simulated tempering */
96 real simtemp_high; /* the high temperature for simulated tempering */
97 real *temperatures; /* the range of temperatures used for simulated tempering */
98 } t_simtemp;
100 typedef struct t_lambda {
101 int nstdhdl; /* The frequency for calculating dhdl */
102 double init_lambda; /* fractional value of lambda (usually will use
103 init_fep_state, this will only be for slow growth,
104 and for legacy free energy code. Only has a
105 valid value if positive) */
106 int init_fep_state; /* the initial number of the state */
107 double delta_lambda; /* change of lambda per time step (fraction of (0.1) */
108 int edHdLPrintEnergy; /* print no, total or potential energies in dhdl */
109 int n_lambda; /* The number of foreign lambda points */
110 double **all_lambda; /* The array of all lambda values */
111 int lambda_neighbors; /* The number of neighboring lambda states to
112 calculate the energy for in up and down directions
113 (-1 for all) */
114 int lambda_start_n; /* The first lambda to calculate energies for */
115 int lambda_stop_n; /* The last lambda +1 to calculate energies for */
116 real sc_alpha; /* free energy soft-core parameter */
117 int sc_power; /* lambda power for soft-core interactions */
118 real sc_r_power; /* r power for soft-core interactions */
119 real sc_sigma; /* free energy soft-core sigma when c6 or c12=0 */
120 real sc_sigma_min; /* free energy soft-core sigma for ????? */
121 gmx_bool bScCoul; /* use softcore for the coulomb portion as well (default FALSE) */
122 gmx_bool separate_dvdl[efptNR]; /* whether to print the dvdl term associated with
123 this term; if it is not specified as separate,
124 it is lumped with the FEP term */
125 int separate_dhdl_file; /* whether to write a separate dhdl.xvg file
126 note: NOT a gmx_bool, but an enum */
127 int dhdl_derivatives; /* whether to calculate+write dhdl derivatives
128 note: NOT a gmx_bool, but an enum */
129 int dh_hist_size; /* The maximum table size for the dH histogram */
130 double dh_hist_spacing; /* The spacing for the dH histogram */
131 } t_lambda;
133 typedef struct t_expanded {
134 int nstexpanded; /* The frequency of expanded ensemble state changes */
135 int elamstats; /* which type of move updating do we use for lambda monte carlo (or no for none) */
136 int elmcmove; /* what move set will be we using for state space moves */
137 int elmceq; /* the method we use to decide of we have equilibrated the weights */
138 int equil_n_at_lam; /* the minumum number of samples at each lambda for deciding whether we have reached a minimum */
139 real equil_wl_delta; /* WL delta at which we stop equilibrating weights */
140 real equil_ratio; /* use the ratio of weights (ratio of minimum to maximum) to decide when to stop equilibrating */
141 int equil_steps; /* after equil_steps steps we stop equilibrating the weights */
142 int equil_samples; /* after equil_samples total samples (steps/nstfep), we stop equilibrating the weights */
143 int lmc_seed; /* random number seed for lambda mc switches */
144 gmx_bool minvar; /* whether to use minumum variance weighting */
145 int minvarmin; /* the number of samples needed before kicking into minvar routine */
146 real minvar_const; /* the offset for the variance in MinVar */
147 int c_range; /* range of cvalues used for BAR */
148 gmx_bool bSymmetrizedTMatrix; /* whether to print symmetrized matrices */
149 int nstTij; /* How frequently to print the transition matrices */
150 int lmc_repeats; /* number of repetitions in the MC lambda jumps */ /*MRS -- VERIFY THIS */
151 int lmc_forced_nstart; /* minimum number of samples for each state before free sampling */ /* MRS -- VERIFY THIS! */
152 int gibbsdeltalam; /* distance in lambda space for the gibbs interval */
153 real wl_scale; /* scaling factor for wang-landau */
154 real wl_ratio; /* ratio between largest and smallest number for freezing the weights */
155 real init_wl_delta; /* starting delta for wang-landau */
156 gmx_bool bWLoneovert; /* use one over t convergece for wang-landau when the delta get sufficiently small */
157 gmx_bool bInit_weights; /* did we initialize the weights? TODO: REMOVE FOR 5.0, no longer needed with new logic */
158 real mc_temp; /* To override the main temperature, or define it if it's not defined */
159 real *init_lambda_weights; /* user-specified initial weights to start with */
160 } t_expanded;
163 /* Abstract types for enforced rotation only defined in pull_rotation.c */
164 typedef struct gmx_enfrot *gmx_enfrot_t;
165 typedef struct gmx_enfrotgrp *gmx_enfrotgrp_t;
167 typedef struct {
168 int eType; /* Rotation type for this group */
169 int bMassW; /* Use mass-weighed positions? */
170 int nat; /* Number of atoms in the group */
171 int *ind; /* The global atoms numbers */
172 rvec *x_ref; /* The reference positions */
173 rvec vec; /* The normalized rotation vector */
174 real rate; /* Rate of rotation (degree/ps) */
175 real k; /* Force constant (kJ/(mol nm^2) */
176 rvec pivot; /* Pivot point of rotation axis (nm) */
177 int eFittype; /* Type of fit to determine actual group angle */
178 int PotAngle_nstep; /* Number of angles around the reference angle
179 for which the rotation potential is also
180 evaluated (for fit type 'potential' only) */
181 real PotAngle_step; /* Distance between two angles in degrees (for
182 fit type 'potential' only) */
183 real slab_dist; /* Slab distance (nm) */
184 real min_gaussian; /* Minimum value the gaussian must have so that
185 the force is actually evaluated */
186 real eps; /* Additive constant for radial motion2 and
187 flexible2 potentials (nm^2) */
188 gmx_enfrotgrp_t enfrotgrp; /* Stores non-inputrec rotation data per group */
189 } t_rotgrp;
191 typedef struct t_rot {
192 int ngrp; /* Number of rotation groups */
193 int nstrout; /* Output frequency for main rotation outfile */
194 int nstsout; /* Output frequency for per-slab data */
195 t_rotgrp *grp; /* Groups to rotate */
196 gmx_enfrot_t enfrot; /* Stores non-inputrec enforced rotation data */
197 } t_rot;
199 /* Abstract type for IMD only defined in IMD.c */
200 struct t_gmx_IMD;
202 typedef struct t_IMD {
203 int nat; /* Number of interactive atoms */
204 int *ind; /* The global indices of the interactive atoms */
205 struct t_gmx_IMD *setup; /* Stores non-inputrec IMD data */
206 } t_IMD;
208 /* Abstract types for position swapping only defined in swapcoords.cpp */
209 typedef struct t_swap *gmx_swapcoords_t;
211 typedef struct t_swapGroup {
212 char *molname; /* Name of the swap group, e.g. NA, CL, SOL */
213 int nat; /* Number of atoms in this group */
214 int *ind; /* The global ion group atoms numbers */
215 int nmolReq[eCompNR]; /* Requested number of molecules of this type
216 per compartment */
217 } t_swapGroup;
219 typedef struct t_swapcoords {
220 int nstswap; /* Every how many steps a swap is attempted? */
221 gmx_bool massw_split[2]; /* Use mass-weighted positions in split group? */
222 real cyl0r, cyl1r; /* Split cylinders defined by radius, upper and */
223 real cyl0u, cyl1u; /* ... lower extension. The split cylinders de- */
224 real cyl0l, cyl1l; /* ... fine the channels and are each anchored */
225 /* ... in the center of the split group */
226 int nAverage; /* Coupling constant (nr of swap attempt steps) */
227 real threshold; /* Ion counts may deviate from the requested
228 values by +-threshold before a swap is done */
229 real bulkOffset[eCompNR]; /* Offset of the swap layer (='bulk') w.r.t.
230 the compartment-defining layers */
231 int ngrp; /* Number of groups to be controlled */
232 t_swapGroup *grp; /* All swap groups, including split and solvent */
233 gmx_swapcoords_t si_priv; /* swap private data accessible in
234 * swapcoords.cpp */
235 } t_swapcoords;
237 struct t_inputrec
239 t_inputrec();
240 explicit t_inputrec(const t_inputrec &) = delete;
241 t_inputrec &operator=(const t_inputrec &) = delete;
242 ~t_inputrec();
244 int eI; /* Integration method */
245 gmx_int64_t nsteps; /* number of steps to be taken */
246 int simulation_part; /* Used in checkpointing to separate chunks */
247 gmx_int64_t init_step; /* start at a stepcount >0 (used w. convert-tpr) */
248 int nstcalcenergy; /* frequency of energy calc. and T/P coupl. upd. */
249 int cutoff_scheme; /* group or verlet cutoffs */
250 int ns_type; /* which ns method should we use? */
251 int nstlist; /* number of steps before pairlist is generated */
252 int ndelta; /* number of cells per rlong */
253 int nstcomm; /* number of steps after which center of mass */
254 /* motion is removed */
255 int comm_mode; /* Center of mass motion removal algorithm */
256 int nstlog; /* number of steps after which print to logfile */
257 int nstxout; /* number of steps after which X is output */
258 int nstvout; /* id. for V */
259 int nstfout; /* id. for F */
260 int nstenergy; /* number of steps after which energies printed */
261 int nstxout_compressed; /* id. for compressed trj (.xtc,.tng) */
262 double init_t; /* initial time (ps) */
263 double delta_t; /* time step (ps) */
264 real x_compression_precision; /* precision of x in compressed trajectory file */
265 real fourier_spacing; /* requested fourier_spacing, when nk? not set */
266 int nkx, nky, nkz; /* number of k vectors in each spatial dimension*/
267 /* for fourier methods for long range electrost.*/
268 int pme_order; /* interpolation order for PME */
269 real ewald_rtol; /* Real space tolerance for Ewald, determines */
270 /* the real/reciprocal space relative weight */
271 real ewald_rtol_lj; /* Real space tolerance for LJ-Ewald */
272 int ewald_geometry; /* normal/3d ewald, or pseudo-2d LR corrections */
273 real epsilon_surface; /* Epsilon for PME dipole correction */
274 int ljpme_combination_rule; /* Type of combination rule in LJ-PME */
275 int ePBC; /* Type of periodic boundary conditions */
276 int bPeriodicMols; /* Periodic molecules */
277 gmx_bool bContinuation; /* Continuation run: starting state is correct */
278 int etc; /* temperature coupling */
279 int nsttcouple; /* interval in steps for temperature coupling */
280 gmx_bool bPrintNHChains; /* whether to print nose-hoover chains */
281 int epc; /* pressure coupling */
282 int epct; /* pressure coupling type */
283 int nstpcouple; /* interval in steps for pressure coupling */
284 real tau_p; /* pressure coupling time (ps) */
285 tensor ref_p; /* reference pressure (kJ/(mol nm^3)) */
286 tensor compress; /* compressability ((mol nm^3)/kJ) */
287 int refcoord_scaling; /* How to scale absolute reference coordinates */
288 rvec posres_com; /* The COM of the posres atoms */
289 rvec posres_comB; /* The B-state COM of the posres atoms */
290 int andersen_seed; /* Random seed for Andersen thermostat (obsolete) */
291 real verletbuf_tol; /* Per atom pair energy drift tolerance (kJ/mol/ps/atom) for list buffer */
292 real rlist; /* short range pairlist cut-off (nm) */
293 real rtpi; /* Radius for test particle insertion */
294 int coulombtype; /* Type of electrostatics treatment */
295 int coulomb_modifier; /* Modify the Coulomb interaction */
296 real rcoulomb_switch; /* Coulomb switch range start (nm) */
297 real rcoulomb; /* Coulomb cutoff (nm) */
298 real epsilon_r; /* relative dielectric constant */
299 real epsilon_rf; /* relative dielectric constant of the RF */
300 bool implicit_solvent; /* Always false (no longer supported */
301 int vdwtype; /* Type of Van der Waals treatment */
302 int vdw_modifier; /* Modify the VdW interaction */
303 real rvdw_switch; /* Van der Waals switch range start (nm) */
304 real rvdw; /* Van der Waals cutoff (nm) */
305 int eDispCorr; /* Perform Long range dispersion corrections */
306 real tabext; /* Extension of the table beyond the cut-off, *
307 * as well as the table length for 1-4 interac. */
308 real shake_tol; /* tolerance for shake */
309 int efep; /* free energy calculations */
310 t_lambda *fepvals; /* Data for the FEP state */
311 gmx_bool bSimTemp; /* Whether to do simulated tempering */
312 t_simtemp *simtempvals; /* Variables for simulated tempering */
313 gmx_bool bExpanded; /* Whether expanded ensembles are used */
314 t_expanded *expandedvals; /* Expanded ensemble parameters */
315 int eDisre; /* Type of distance restraining */
316 real dr_fc; /* force constant for ta_disre */
317 int eDisreWeighting; /* type of weighting of pairs in one restraints */
318 gmx_bool bDisreMixed; /* Use comb of time averaged and instan. viol's */
319 int nstdisreout; /* frequency of writing pair distances to enx */
320 real dr_tau; /* time constant for memory function in disres */
321 real orires_fc; /* force constant for orientational restraints */
322 real orires_tau; /* time constant for memory function in orires */
323 int nstorireout; /* frequency of writing tr(SD) to enx */
324 real em_stepsize; /* The stepsize for updating */
325 real em_tol; /* The tolerance */
326 int niter; /* Number of iterations for convergence of */
327 /* steepest descent in relax_shells */
328 real fc_stepsize; /* Stepsize for directional minimization */
329 /* in relax_shells */
330 int nstcgsteep; /* number of steps after which a steepest */
331 /* descents step is done while doing cg */
332 int nbfgscorr; /* Number of corrections to the hessian to keep */
333 int eConstrAlg; /* Type of constraint algorithm */
334 int nProjOrder; /* Order of the LINCS Projection Algorithm */
335 real LincsWarnAngle; /* If bond rotates more than %g degrees, warn */
336 int nLincsIter; /* Number of iterations in the final Lincs step */
337 gmx_bool bShakeSOR; /* Use successive overrelaxation for shake */
338 real bd_fric; /* Friction coefficient for BD (amu/ps) */
339 gmx_int64_t ld_seed; /* Random seed for SD and BD */
340 int nwall; /* The number of walls */
341 int wall_type; /* The type of walls */
342 real wall_r_linpot; /* The potentail is linear for r<=wall_r_linpot */
343 int wall_atomtype[2]; /* The atom type for walls */
344 real wall_density[2]; /* Number density for walls */
345 real wall_ewald_zfac; /* Scaling factor for the box for Ewald */
347 /* COM pulling data */
348 gmx_bool bPull; /* Do we do COM pulling? */
349 struct pull_params_t *pull; /* The data for center of mass pulling */
350 // TODO: Remove this by converting pull into a ForceProvider
351 struct pull_t *pull_work; /* The COM pull force calculation data structure */
353 /* AWH bias data */
354 gmx_bool bDoAwh; /* Use awh biasing for PMF calculations? */
355 gmx::AwhParams *awhParams; /* AWH biasing parameters */
357 /* Enforced rotation data */
358 gmx_bool bRot; /* Calculate enforced rotation potential(s)? */
359 t_rot *rot; /* The data for enforced rotation potentials */
361 int eSwapCoords; /* Do ion/water position exchanges (CompEL)? */
362 t_swapcoords *swap;
364 gmx_bool bIMD; /* Allow interactive MD sessions for this .tpr? */
365 t_IMD *imd; /* Interactive molecular dynamics */
367 real cos_accel; /* Acceleration for viscosity calculation */
368 tensor deform; /* Triclinic deformation velocities (nm/ps) */
369 int userint1; /* User determined parameters */
370 int userint2;
371 int userint3;
372 int userint4;
373 real userreal1;
374 real userreal2;
375 real userreal3;
376 real userreal4;
377 t_grpopts opts; /* Group options */
378 gmx_bool bQMMM; /* QM/MM calculation */
379 int QMconstraints; /* constraints on QM bonds */
380 int QMMMscheme; /* Scheme: ONIOM or normal */
381 real scalefactor; /* factor for scaling the MM charges in QM calc.*/
383 /* Fields for removed features go here (better caching) */
384 gmx_bool bAdress; // Whether AdResS is enabled - always false if a valid .tpr was read
385 gmx_bool useTwinRange; // Whether twin-range scheme is active - always false if a valid .tpr was read
387 gmx::KeyValueTreeObject *params;
390 int ir_optimal_nstcalcenergy(const t_inputrec *ir);
392 int tcouple_min_integration_steps(int etc);
394 int ir_optimal_nsttcouple(const t_inputrec *ir);
396 int pcouple_min_integration_steps(int epc);
398 int ir_optimal_nstpcouple(const t_inputrec *ir);
400 /* Returns if the Coulomb force or potential is switched to zero */
401 gmx_bool ir_coulomb_switched(const t_inputrec *ir);
403 /* Returns if the Coulomb interactions are zero beyond the rcoulomb.
404 * Note: always returns TRUE for the Verlet cut-off scheme.
406 gmx_bool ir_coulomb_is_zero_at_cutoff(const t_inputrec *ir);
408 /* As ir_coulomb_is_zero_at_cutoff, but also returns TRUE for user tabulated
409 * interactions, since these might be zero beyond rcoulomb.
411 gmx_bool ir_coulomb_might_be_zero_at_cutoff(const t_inputrec *ir);
413 /* Returns if the Van der Waals force or potential is switched to zero */
414 gmx_bool ir_vdw_switched(const t_inputrec *ir);
416 /* Returns if the Van der Waals interactions are zero beyond the rvdw.
417 * Note: always returns TRUE for the Verlet cut-off scheme.
419 gmx_bool ir_vdw_is_zero_at_cutoff(const t_inputrec *ir);
421 /* As ir_vdw_is_zero_at_cutoff, but also returns TRUE for user tabulated
422 * interactions, since these might be zero beyond rvdw.
424 gmx_bool ir_vdw_might_be_zero_at_cutoff(const t_inputrec *ir);
426 /*! \brief Free memory from input record.
428 * All arrays and pointers will be freed.
430 * \param[in] ir The data structure
432 void done_inputrec(t_inputrec *ir);
434 void pr_inputrec(FILE *fp, int indent, const char *title, const t_inputrec *ir,
435 gmx_bool bMDPformat);
437 void cmp_inputrec(FILE *fp, const t_inputrec *ir1, const t_inputrec *ir2, real ftol, real abstol);
439 void comp_pull_AB(FILE *fp, pull_params_t *pull, real ftol, real abstol);
442 gmx_bool inputrecDeform(const t_inputrec *ir);
444 gmx_bool inputrecDynamicBox(const t_inputrec *ir);
446 gmx_bool inputrecPreserveShape(const t_inputrec *ir);
448 gmx_bool inputrecNeedMutot(const t_inputrec *ir);
450 gmx_bool inputrecTwinRange(const t_inputrec *ir);
452 gmx_bool inputrecExclForces(const t_inputrec *ir);
454 gmx_bool inputrecNptTrotter(const t_inputrec *ir);
456 gmx_bool inputrecNvtTrotter(const t_inputrec *ir);
458 gmx_bool inputrecNphTrotter(const t_inputrec *ir);
460 /*! \brief Return true if the simulation is 2D periodic with two walls. */
461 bool inputrecPbcXY2Walls(const t_inputrec *ir);
463 /*! \brief Returns true for MD integator with T and/or P-coupling that supports
464 * calculating the conserved energy quantity.
466 bool integratorHasConservedEnergyQuantity(const t_inputrec *ir);
468 /*! \brief Returns true when temperature is coupled or constant. */
469 bool integratorHasReferenceTemperature(const t_inputrec *ir);
471 /*! \brief Return the number of bounded dimensions
473 * \param[in] ir The input record with MD parameters
474 * \return the number of dimensions in which
475 * the coordinates of the particles are bounded, starting at X.
477 int inputrec2nboundeddim(const t_inputrec *ir);
479 /*! \brief Returns the number of degrees of freedom in center of mass motion
481 * \param[in] ir the inputrec structure
482 * \return the number of degrees of freedom of the center of mass
484 int ndof_com(const t_inputrec *ir);
486 #endif /* GMX_MDTYPES_INPUTREC_H */