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 * check out http://www.gromacs.org for more information.
7 * Copyright (c) 2012,2013, by the GROMACS development team, led by
8 * David van der Spoel, Berk Hess, Erik Lindahl, and including many
9 * others, as listed in the AUTHORS file in the top-level source
10 * 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.
43 #include "../sysstuff.h"
51 int n
; /* Number of terms */
52 real
*a
; /* Coeffients (V / nm ) */
53 real
*phi
; /* Phase angles */
57 real E0
; /* Field strength (V/nm) */
58 real omega
; /* Frequency (1/ps) */
59 real t0
; /* Centre of the Gaussian pulse (ps) */
60 real sigma
; /* Width of the Gaussian pulse (FWHM) (ps) */
63 #define EGP_EXCL (1<<0)
64 #define EGP_TABLE (1<<1)
67 int ngtc
; /* # T-Coupl groups */
68 int nhchainlength
; /* # of nose-hoover chains per group */
69 int ngacc
; /* # Accelerate groups */
70 int ngfrz
; /* # Freeze groups */
71 int ngener
; /* # Ener groups */
72 real
*nrdf
; /* Nr of degrees of freedom in a group */
73 real
*ref_t
; /* Coupling temperature per group */
74 int *annealing
; /* No/simple/periodic SA for each group */
75 int *anneal_npoints
; /* Number of annealing time points per grp */
76 real
**anneal_time
; /* For ea. group: Time points */
77 real
**anneal_temp
; /* For ea. grp: Temperature at these times */
78 /* Final temp after all intervals is ref_t */
79 real
*tau_t
; /* Tau coupling time */
80 rvec
*acc
; /* Acceleration per group */
81 ivec
*nFreeze
; /* Freeze the group in each direction ? */
82 int *egp_flags
; /* Exclusions/tables of energy group pairs */
85 int ngQM
; /* nr of QM groups */
86 int *QMmethod
; /* Level of theory in the QM calculation */
87 int *QMbasis
; /* Basisset in the QM calculation */
88 int *QMcharge
; /* Total charge in the QM region */
89 int *QMmult
; /* Spin multiplicicty in the QM region */
90 gmx_bool
*bSH
; /* surface hopping (diabatic hop only) */
91 int *CASorbitals
; /* number of orbiatls in the active space */
92 int *CASelectrons
; /* number of electrons in the active space */
93 real
*SAon
; /* at which gap (A.U.) the SA is switched on */
95 int *SAsteps
; /* in how many steps SA goes from 1-1 to 0.5-0.5*/
101 epgrppbcNONE
, epgrppbcREFAT
, epgrppbcCOS
105 int nat
; /* Number of atoms in the pull group */
106 atom_id
*ind
; /* The global atoms numbers */
107 int nat_loc
; /* Number of local pull atoms */
108 int nalloc_loc
; /* Allocation size for ind_loc and weight_loc */
109 atom_id
*ind_loc
; /* Local pull indices */
110 int nweight
; /* The number of weights (0 or nat) */
111 real
*weight
; /* Weights (use all 1 when weight==NULL) */
112 real
*weight_loc
; /* Weights for the local indices */
113 int epgrppbc
; /* The type of pbc for this pull group, see enum above */
114 atom_id pbcatom
; /* The reference atom for pbc (global number) */
115 rvec vec
; /* The pull vector, direction or position */
116 rvec init
; /* Initial reference displacement */
117 real rate
; /* Rate of motion (nm/ps) */
118 real k
; /* force constant */
119 real kB
; /* force constant for state B */
120 real wscale
; /* scaling factor for the weights: sum w m/sum w w m */
121 real invtm
; /* inverse total mass of the group: 1/wscale sum w m */
122 dvec x
; /* center of mass before update */
123 dvec xp
; /* center of mass after update before constraining */
124 dvec dr
; /* The distance from the reference group */
125 double f_scal
; /* Scalar force for directional pulling */
126 dvec f
; /* force due to the pulling/constraining */
130 int eSimTempScale
; /* simulated temperature scaling; linear or exponential */
131 real simtemp_low
; /* the low temperature for simulated tempering */
132 real simtemp_high
; /* the high temperature for simulated tempering */
133 real
*temperatures
; /* the range of temperatures used for simulated tempering */
137 int nstdhdl
; /* The frequency for calculating dhdl */
138 double init_lambda
; /* fractional value of lambda (usually will use
139 init_fep_state, this will only be for slow growth,
140 and for legacy free energy code. Only has a
141 valid value if positive) */
142 int init_fep_state
; /* the initial number of the state */
143 double delta_lambda
; /* change of lambda per time step (fraction of (0.1) */
144 gmx_bool bPrintEnergy
; /* Whether to print the energy in the dhdl */
145 int n_lambda
; /* The number of foreign lambda points */
146 double **all_lambda
; /* The array of all lambda values */
147 int lambda_neighbors
; /* The number of neighboring lambda states to
148 calculate the energy for in up and down directions
150 int lambda_start_n
; /* The first lambda to calculate energies for */
151 int lambda_stop_n
; /* The last lambda +1 to calculate energies for */
152 real sc_alpha
; /* free energy soft-core parameter */
153 int sc_power
; /* lambda power for soft-core interactions */
154 real sc_r_power
; /* r power for soft-core interactions */
155 real sc_sigma
; /* free energy soft-core sigma when c6 or c12=0 */
156 real sc_sigma_min
; /* free energy soft-core sigma for ????? */
157 gmx_bool bScCoul
; /* use softcore for the coulomb portion as well (default FALSE) */
158 gmx_bool separate_dvdl
[efptNR
]; /* whether to print the dvdl term associated with
159 this term; if it is not specified as separate,
160 it is lumped with the FEP term */
161 int separate_dhdl_file
; /* whether to write a separate dhdl.xvg file
162 note: NOT a gmx_bool, but an enum */
163 int dhdl_derivatives
; /* whether to calculate+write dhdl derivatives
164 note: NOT a gmx_bool, but an enum */
165 int dh_hist_size
; /* The maximum table size for the dH histogram */
166 double dh_hist_spacing
; /* The spacing for the dH histogram */
170 int nstexpanded
; /* The frequency of expanded ensemble state changes */
171 int elamstats
; /* which type of move updating do we use for lambda monte carlo (or no for none) */
172 int elmcmove
; /* what move set will be we using for state space moves */
173 int elmceq
; /* the method we use to decide of we have equilibrated the weights */
174 int equil_n_at_lam
; /* the minumum number of samples at each lambda for deciding whether we have reached a minimum */
175 real equil_wl_delta
; /* WL delta at which we stop equilibrating weights */
176 real equil_ratio
; /* use the ratio of weights (ratio of minimum to maximum) to decide when to stop equilibrating */
177 int equil_steps
; /* after equil_steps steps we stop equilibrating the weights */
178 int equil_samples
; /* after equil_samples total samples (steps/nstfep), we stop equilibrating the weights */
179 int lmc_seed
; /* random number seed for lambda mc switches */
180 gmx_bool minvar
; /* whether to use minumum variance weighting */
181 int minvarmin
; /* the number of samples needed before kicking into minvar routine */
182 real minvar_const
; /* the offset for the variance in MinVar */
183 int c_range
; /* range of cvalues used for BAR */
184 gmx_bool bSymmetrizedTMatrix
; /* whether to print symmetrized matrices */
185 int nstTij
; /* How frequently to print the transition matrices */
186 int lmc_repeats
; /* number of repetitions in the MC lambda jumps */ /*MRS -- VERIFY THIS */
187 int lmc_forced_nstart
; /* minimum number of samples for each state before free sampling */ /* MRS -- VERIFY THIS! */
188 int gibbsdeltalam
; /* distance in lambda space for the gibbs interval */
189 real wl_scale
; /* scaling factor for wang-landau */
190 real wl_ratio
; /* ratio between largest and smallest number for freezing the weights */
191 real init_wl_delta
; /* starting delta for wang-landau */
192 gmx_bool bWLoneovert
; /* use one over t convergece for wang-landau when the delta get sufficiently small */
193 gmx_bool bInit_weights
; /* did we initialize the weights? */
194 real mc_temp
; /* To override the main temperature, or define it if it's not defined */
195 real
*init_lambda_weights
; /* user-specified initial weights to start with */
199 int ngrp
; /* number of groups */
200 int eGeom
; /* pull geometry */
201 ivec dim
; /* used to select components for constraint */
202 real cyl_r1
; /* radius of cylinder for dynamic COM */
203 real cyl_r0
; /* radius of cylinder including switch length */
204 real constr_tol
; /* absolute tolerance for constraints in (nm) */
205 int nstxout
; /* Output frequency for pull x */
206 int nstfout
; /* Output frequency for pull f */
207 int ePBC
; /* the boundary conditions */
208 int npbcdim
; /* do pbc in dims 0 <= dim < npbcdim */
209 gmx_bool bRefAt
; /* do we need reference atoms for a group COM ? */
210 int cosdim
; /* dimension for cosine weighting, -1 if none */
211 gmx_bool bVirial
; /* do we need to add the pull virial? */
212 t_pullgrp
*grp
; /* groups to pull/restrain/etc/ */
213 t_pullgrp
*dyna
; /* dynamic groups for use with local constraints */
214 rvec
*rbuf
; /* COM calculation buffer */
215 dvec
*dbuf
; /* COM calculation buffer */
216 double *dbuf_cyl
; /* cylinder ref. groups COM calculation buffer */
218 FILE *out_x
; /* output file for pull data */
219 FILE *out_f
; /* output file for pull data */
223 /* Abstract types for enforced rotation only defined in pull_rotation.c */
224 typedef struct gmx_enfrot
*gmx_enfrot_t
;
225 typedef struct gmx_enfrotgrp
*gmx_enfrotgrp_t
;
228 int eType
; /* Rotation type for this group */
229 int bMassW
; /* Use mass-weighed positions? */
230 int nat
; /* Number of atoms in the group */
231 atom_id
*ind
; /* The global atoms numbers */
232 rvec
*x_ref
; /* The reference positions */
233 rvec vec
; /* The normalized rotation vector */
234 real rate
; /* Rate of rotation (degree/ps) */
235 real k
; /* Force constant (kJ/(mol nm^2) */
236 rvec pivot
; /* Pivot point of rotation axis (nm) */
237 int eFittype
; /* Type of fit to determine actual group angle */
238 int PotAngle_nstep
; /* Number of angles around the reference angle
239 for which the rotation potential is also
240 evaluated (for fit type 'potential' only) */
241 real PotAngle_step
; /* Distance between two angles in degrees (for
242 fit type 'potential' only) */
243 real slab_dist
; /* Slab distance (nm) */
244 real min_gaussian
; /* Minimum value the gaussian must have so that
245 the force is actually evaluated */
246 real eps
; /* Additive constant for radial motion2 and
247 flexible2 potentials (nm^2) */
248 gmx_enfrotgrp_t enfrotgrp
; /* Stores non-inputrec rotation data per group */
252 int ngrp
; /* Number of rotation groups */
253 int nstrout
; /* Output frequency for main rotation outfile */
254 int nstsout
; /* Output frequency for per-slab data */
255 t_rotgrp
*grp
; /* Groups to rotate */
256 gmx_enfrot_t enfrot
; /* Stores non-inputrec enforced rotation data */
261 int type
; /* type of AdResS simulation */
262 gmx_bool bnew_wf
; /* enable new AdResS weighting function */
263 gmx_bool bchempot_dx
; /*true:interaction table format input is F=-dmu/dx false: dmu_dwp */
264 gmx_bool btf_full_box
; /* true: appy therm force everywhere in the box according to table false: only in hybrid region */
265 real const_wf
; /* value of weighting function for eAdressConst */
266 real ex_width
; /* center of the explicit zone */
267 real hy_width
; /* width of the hybrid zone */
268 int icor
; /* type of interface correction */
269 int site
; /* AdResS CG site location */
270 rvec refs
; /* Coordinates for AdResS reference */
271 real ex_forcecap
; /* in the hybrid zone, cap forces large then this to adress_ex_forcecap */
272 gmx_bool do_hybridpairs
; /* If true pair interaction forces are also scaled in an adress way*/
274 int * tf_table_index
; /* contains mapping of energy group index -> i-th adress tf table*/
281 int eI
; /* Integration method */
282 gmx_large_int_t nsteps
; /* number of steps to be taken */
283 int simulation_part
; /* Used in checkpointing to separate chunks */
284 gmx_large_int_t init_step
; /* start at a stepcount >0 (used w. tpbconv) */
285 int nstcalcenergy
; /* frequency of energy calc. and T/P coupl. upd. */
286 int cutoff_scheme
; /* group or verlet cutoffs */
287 int ns_type
; /* which ns method should we use? */
288 int nstlist
; /* number of steps before pairlist is generated */
289 int ndelta
; /* number of cells per rlong */
290 int nstcomm
; /* number of steps after which center of mass */
291 /* motion is removed */
292 int comm_mode
; /* Center of mass motion removal algorithm */
293 int nstcheckpoint
; /* checkpointing frequency */
294 int nstlog
; /* number of steps after which print to logfile */
295 int nstxout
; /* number of steps after which X is output */
296 int nstvout
; /* id. for V */
297 int nstfout
; /* id. for F */
298 int nstenergy
; /* number of steps after which energies printed */
299 int nstxtcout
; /* id. for compressed trj (.xtc) */
300 double init_t
; /* initial time (ps) */
301 double delta_t
; /* time step (ps) */
302 real xtcprec
; /* precision of xtc file */
303 real fourier_spacing
; /* requested fourier_spacing, when nk? not set */
304 int nkx
, nky
, nkz
; /* number of k vectors in each spatial dimension*/
305 /* for fourier methods for long range electrost.*/
306 int pme_order
; /* interpolation order for PME */
307 real ewald_rtol
; /* Real space tolerance for Ewald, determines */
308 /* the real/reciprocal space relative weight */
309 int ewald_geometry
; /* normal/3d ewald, or pseudo-2d LR corrections */
310 real epsilon_surface
; /* Epsilon for PME dipole correction */
311 gmx_bool bOptFFT
; /* optimize the fft plan at start */
312 int ePBC
; /* Type of periodic boundary conditions */
313 int bPeriodicMols
; /* Periodic molecules */
314 gmx_bool bContinuation
; /* Continuation run: starting state is correct */
315 int etc
; /* temperature coupling */
316 int nsttcouple
; /* interval in steps for temperature coupling */
317 gmx_bool bPrintNHChains
; /* whether to print nose-hoover chains */
318 int epc
; /* pressure coupling */
319 int epct
; /* pressure coupling type */
320 int nstpcouple
; /* interval in steps for pressure coupling */
321 real tau_p
; /* pressure coupling time (ps) */
322 tensor ref_p
; /* reference pressure (kJ/(mol nm^3)) */
323 tensor compress
; /* compressability ((mol nm^3)/kJ) */
324 int refcoord_scaling
; /* How to scale absolute reference coordinates */
325 rvec posres_com
; /* The COM of the posres atoms */
326 rvec posres_comB
; /* The B-state COM of the posres atoms */
327 int andersen_seed
; /* Random seed for Andersen thermostat (obsolete) */
328 real verletbuf_drift
; /* Max. drift (kJ/mol/ps/atom) for list buffer */
329 real rlist
; /* short range pairlist cut-off (nm) */
330 real rlistlong
; /* long range pairlist cut-off (nm) */
331 int nstcalclr
; /* Frequency of evaluating direct space long-range interactions */
332 real rtpi
; /* Radius for test particle insertion */
333 int coulombtype
; /* Type of electrostatics treatment */
334 int coulomb_modifier
; /* Modify the Coulomb interaction */
335 real rcoulomb_switch
; /* Coulomb switch range start (nm) */
336 real rcoulomb
; /* Coulomb cutoff (nm) */
337 real epsilon_r
; /* relative dielectric constant */
338 real epsilon_rf
; /* relative dielectric constant of the RF */
339 int implicit_solvent
; /* No (=explicit water), or GBSA solvent models */
340 int gb_algorithm
; /* Algorithm to use for calculation Born radii */
341 int nstgbradii
; /* Frequency of updating Generalized Born radii */
342 real rgbradii
; /* Cutoff for GB radii calculation */
343 real gb_saltconc
; /* Salt concentration (M) for GBSA models */
344 real gb_epsilon_solvent
; /* dielectric coeff. of implicit solvent */
345 real gb_obc_alpha
; /* 1st scaling factor for Bashford-Case GB */
346 real gb_obc_beta
; /* 2nd scaling factor for Bashford-Case GB */
347 real gb_obc_gamma
; /* 3rd scaling factor for Bashford-Case GB */
348 real gb_dielectric_offset
; /* Dielectric offset for Still/HCT/OBC */
349 int sa_algorithm
; /* Algorithm for SA part of GBSA */
350 real sa_surface_tension
; /* Energy factor for SA part of GBSA */
351 int vdwtype
; /* Type of Van der Waals treatment */
352 int vdw_modifier
; /* Modify the VdW interaction */
353 real rvdw_switch
; /* Van der Waals switch range start (nm) */
354 real rvdw
; /* Van der Waals cutoff (nm) */
355 int eDispCorr
; /* Perform Long range dispersion corrections */
356 real tabext
; /* Extension of the table beyond the cut-off, *
357 * as well as the table length for 1-4 interac. */
358 real shake_tol
; /* tolerance for shake */
359 int efep
; /* free energy calculations */
360 t_lambda
*fepvals
; /* Data for the FEP state */
361 gmx_bool bSimTemp
; /* Whether to do simulated tempering */
362 t_simtemp
*simtempvals
; /* Variables for simulated tempering */
363 gmx_bool bExpanded
; /* Whether expanded ensembles are used */
364 t_expanded
*expandedvals
; /* Expanded ensemble parameters */
365 int eDisre
; /* Type of distance restraining */
366 real dr_fc
; /* force constant for ta_disre */
367 int eDisreWeighting
; /* type of weighting of pairs in one restraints */
368 gmx_bool bDisreMixed
; /* Use comb of time averaged and instan. viol's */
369 int nstdisreout
; /* frequency of writing pair distances to enx */
370 real dr_tau
; /* time constant for memory function in disres */
371 real orires_fc
; /* force constant for orientational restraints */
372 real orires_tau
; /* time constant for memory function in orires */
373 int nstorireout
; /* frequency of writing tr(SD) to enx */
374 real dihre_fc
; /* force constant for dihedral restraints (obsolete) */
375 real em_stepsize
; /* The stepsize for updating */
376 real em_tol
; /* The tolerance */
377 int niter
; /* Number of iterations for convergence of */
378 /* steepest descent in relax_shells */
379 real fc_stepsize
; /* Stepsize for directional minimization */
380 /* in relax_shells */
381 int nstcgsteep
; /* number of steps after which a steepest */
382 /* descents step is done while doing cg */
383 int nbfgscorr
; /* Number of corrections to the hessian to keep */
384 int eConstrAlg
; /* Type of constraint algorithm */
385 int nProjOrder
; /* Order of the LINCS Projection Algorithm */
386 real LincsWarnAngle
; /* If bond rotates more than %g degrees, warn */
387 int nLincsIter
; /* Number of iterations in the final Lincs step */
388 gmx_bool bShakeSOR
; /* Use successive overrelaxation for shake */
389 real bd_fric
; /* Friction coefficient for BD (amu/ps) */
390 int ld_seed
; /* Random seed for SD and BD */
391 int nwall
; /* The number of walls */
392 int wall_type
; /* The type of walls */
393 real wall_r_linpot
; /* The potentail is linear for r<=wall_r_linpot */
394 int wall_atomtype
[2]; /* The atom type for walls */
395 real wall_density
[2]; /* Number density for walls */
396 real wall_ewald_zfac
; /* Scaling factor for the box for Ewald */
397 int ePull
; /* Type of pulling: no, umbrella or constraint */
398 t_pull
*pull
; /* The data for center of mass pulling */
399 gmx_bool bRot
; /* Calculate enforced rotation potential(s)? */
400 t_rot
*rot
; /* The data for enforced rotation potentials */
401 real cos_accel
; /* Acceleration for viscosity calculation */
402 tensor deform
; /* Triclinic deformation velocities (nm/ps) */
403 int userint1
; /* User determined parameters */
411 t_grpopts opts
; /* Group options */
412 t_cosines ex
[DIM
]; /* Electric field stuff (spatial part) */
413 t_cosines et
[DIM
]; /* Electric field stuff (time part) */
414 gmx_bool bQMMM
; /* QM/MM calculation */
415 int QMconstraints
; /* constraints on QM bonds */
416 int QMMMscheme
; /* Scheme: ONIOM or normal */
417 real scalefactor
; /* factor for scaling the MM charges in QM calc.*/
418 /* parameter needed for AdResS simulation */
419 gmx_bool bAdress
; /* Is AdResS enabled ? */
420 t_adress
*adress
; /* The data for adress simulations */
423 #define DEFORM(ir) ((ir).deform[XX][XX] != 0 || (ir).deform[YY][YY] != 0 || (ir).deform[ZZ][ZZ] != 0 || (ir).deform[YY][XX] != 0 || (ir).deform[ZZ][XX] != 0 || (ir).deform[ZZ][YY] != 0)
425 #define DYNAMIC_BOX(ir) ((ir).epc != epcNO || (ir).eI == eiTPI || DEFORM(ir))
427 #define PRESERVE_SHAPE(ir) ((ir).epc != epcNO && (ir).deform[XX][XX] == 0 && ((ir).epct == epctISOTROPIC || (ir).epct == epctSEMIISOTROPIC))
429 #define NEED_MUTOT(ir) (((ir).coulombtype == eelEWALD || EEL_PME((ir).coulombtype)) && ((ir).ewald_geometry == eewg3DC || (ir).epsilon_surface != 0))
431 #define IR_TWINRANGE(ir) ((ir).rlist > 0 && ((ir).rlistlong == 0 || (ir).rlistlong > (ir).rlist))
433 #define IR_ELEC_FIELD(ir) ((ir).ex[XX].n > 0 || (ir).ex[YY].n > 0 || (ir).ex[ZZ].n > 0)
435 #define IR_EXCL_FORCES(ir) (EEL_FULL((ir).coulombtype) || (EEL_RF((ir).coulombtype) && (ir).coulombtype != eelRF_NEC) || (ir).implicit_solvent != eisNO)
436 /* use pointer definitions of ir here, since that's what's usually used in the code */
437 #define IR_NPT_TROTTER(ir) ((((ir)->eI == eiVV) || ((ir)->eI == eiVVAK)) && (((ir)->epc == epcMTTK) && ((ir)->etc == etcNOSEHOOVER)))
439 #define IR_NVT_TROTTER(ir) ((((ir)->eI == eiVV) || ((ir)->eI == eiVVAK)) && ((!((ir)->epc == epcMTTK)) && ((ir)->etc == etcNOSEHOOVER)))
441 #define IR_NPH_TROTTER(ir) ((((ir)->eI == eiVV) || ((ir)->eI == eiVVAK)) && (((ir)->epc == epcMTTK) && (!(((ir)->etc == etcNOSEHOOVER)))))
