3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
10 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
11 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
12 * Copyright (c) 2001-2004, The GROMACS development team,
13 * check out http://www.gromacs.org for more information.
15 * This program is free software; you can redistribute it and/or
16 * modify it under the terms of the GNU General Public License
17 * as published by the Free Software Foundation; either version 2
18 * of the License, or (at your option) any later version.
20 * If you want to redistribute modifications, please consider that
21 * scientific software is very special. Version control is crucial -
22 * bugs must be traceable. We will be happy to consider code for
23 * inclusion in the official distribution, but derived work must not
24 * be called official GROMACS. Details are found in the README & COPYING
25 * files - if they are missing, get the official version at www.gromacs.org.
27 * To help us fund GROMACS development, we humbly ask that you cite
28 * the papers on the package - you can find them in the top README file.
30 * For more info, check our website at http://www.gromacs.org
33 * GRoups of Organic Molecules in ACtion for Science
45 #include "../sysstuff.h"
53 int n
; /* Number of terms */
54 real
*a
; /* Coeffients (V / nm ) */
55 real
*phi
; /* Phase angles */
59 real E0
; /* Field strength (V/nm) */
60 real omega
; /* Frequency (1/ps) */
61 real t0
; /* Centre of the Gaussian pulse (ps) */
62 real sigma
; /* Width of the Gaussian pulse (FWHM) (ps) */
65 #define EGP_EXCL (1<<0)
66 #define EGP_TABLE (1<<1)
69 int ngtc
; /* # T-Coupl groups */
70 int nhchainlength
; /* # of nose-hoover chains per group */
71 int ngacc
; /* # Accelerate groups */
72 int ngfrz
; /* # Freeze groups */
73 int ngener
; /* # Ener groups */
74 real
*nrdf
; /* Nr of degrees of freedom in a group */
75 real
*ref_t
; /* Coupling temperature per group */
76 int *annealing
; /* No/simple/periodic SA for each group */
77 int *anneal_npoints
; /* Number of annealing time points per grp */
78 real
**anneal_time
; /* For ea. group: Time points */
79 real
**anneal_temp
; /* For ea. grp: Temperature at these times */
80 /* Final temp after all intervals is ref_t */
81 real
*tau_t
; /* Tau coupling time */
82 rvec
*acc
; /* Acceleration per group */
83 ivec
*nFreeze
; /* Freeze the group in each direction ? */
84 int *egp_flags
; /* Exclusions/tables of energy group pairs */
87 int ngQM
; /* nr of QM groups */
88 int *QMmethod
; /* Level of theory in the QM calculation */
89 int *QMbasis
; /* Basisset in the QM calculation */
90 int *QMcharge
; /* Total charge in the QM region */
91 int *QMmult
; /* Spin multiplicicty in the QM region */
92 bool *bSH
; /* surface hopping (diabatic hop only) */
93 int *CASorbitals
; /* number of orbiatls in the active space */
94 int *CASelectrons
;/* number of electrons in the active space */
95 real
*SAon
; /* at which gap (A.U.) the SA is switched on */
97 int *SAsteps
; /* in how many steps SA goes from 1-1 to 0.5-0.5*/
102 enum { 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 ngrp
; /* number of groups */
131 int eGeom
; /* pull geometry */
132 ivec dim
; /* used to select components for constraint */
133 real cyl_r1
; /* radius of cylinder for dynamic COM */
134 real cyl_r0
; /* radius of cylinder including switch length */
135 real constr_tol
; /* absolute tolerance for constraints in (nm) */
136 int nstxout
; /* Output frequency for pull x */
137 int nstfout
; /* Output frequency for pull f */
138 int ePBC
; /* the boundary conditions */
139 int npbcdim
; /* do pbc in dims 0 <= dim < npbcdim */
140 bool bRefAt
; /* do we need reference atoms for a group COM ? */
141 int cosdim
; /* dimension for cosine weighting, -1 if none */
142 bool bVirial
; /* do we need to add the pull virial? */
143 t_pullgrp
*grp
; /* groups to pull/restrain/etc/ */
144 t_pullgrp
*dyna
; /* dynamic groups for use with local constraints */
145 rvec
*rbuf
; /* COM calculation buffer */
146 dvec
*dbuf
; /* COM calculation buffer */
147 double *dbuf_cyl
; /* cylinder ref. groups COM calculation buffer */
149 FILE *out_x
; /* output file for pull data */
150 FILE *out_f
; /* output file for pull data */
154 int eI
; /* Integration method */
155 gmx_large_int_t nsteps
; /* number of steps to be taken */
156 int simulation_part
; /* Used in checkpointing to separate chunks */
157 gmx_large_int_t init_step
; /* start at a stepcount >0 (used w. tpbconv) */
158 int nstcalcenergy
; /* fequency of energy calc. and T/P coupl. upd. */
159 int ns_type
; /* which ns method should we use? */
160 int nstlist
; /* number of steps before pairlist is generated */
161 int ndelta
; /* number of cells per rlong */
162 int nstcomm
; /* number of steps after which center of mass */
163 /* motion is removed */
164 int comm_mode
; /* Center of mass motion removal algorithm */
165 int nstcheckpoint
; /* checkpointing frequency */
166 int nstlog
; /* number of steps after which print to logfile */
167 int nstxout
; /* number of steps after which X is output */
168 int nstvout
; /* id. for V */
169 int nstfout
; /* id. for F */
170 int nstenergy
; /* number of steps after which energies printed */
171 int nstxtcout
; /* id. for compressed trj (.xtc) */
172 double init_t
; /* initial time (ps) */
173 double delta_t
; /* time step (ps) */
174 real xtcprec
; /* precision of xtc file */
175 int nkx
,nky
,nkz
; /* number of k vectors in each spatial dimension*/
176 /* for fourier methods for long range electrost.*/
177 int pme_order
; /* interpolation order for PME */
178 real ewald_rtol
; /* Real space tolerance for Ewald, determines */
179 /* the real/reciprocal space relative weight */
180 int ewald_geometry
; /* normal/3d ewald, or pseudo-2d LR corrections */
181 real epsilon_surface
; /* Epsilon for PME dipole correction */
182 bool bOptFFT
; /* optimize the fft plan at start */
183 int ePBC
; /* Type of periodic boundary conditions */
184 int bPeriodicMols
; /* Periodic molecules */
185 bool bContinuation
; /* Continuation run: starting state is correct */
186 int etc
; /* temperature coupling */
187 int nsttcouple
; /* interval in steps for temperature coupling */
188 int epc
; /* pressure coupling */
189 int epct
; /* pressure coupling type */
190 int nstpcouple
; /* interval in steps for pressure coupling */
191 real tau_p
; /* pressure coupling time (ps) */
192 tensor ref_p
; /* reference pressure (kJ/(mol nm^3)) */
193 tensor compress
; /* compressability ((mol nm^3)/kJ) */
194 int refcoord_scaling
;/* How to scale absolute reference coordinates */
195 rvec posres_com
; /* The COM of the posres atoms */
196 rvec posres_comB
; /* The B-state COM of the posres atoms */
197 int andersen_seed
; /* Random seed for Andersen thermostat. */
198 real rlist
; /* short range pairlist cut-off (nm) */
199 real rlistlong
; /* long range pairlist cut-off (nm) */
200 real rtpi
; /* Radius for test particle insertion */
201 int coulombtype
; /* Type of electrostatics treatment */
202 real rcoulomb_switch
; /* Coulomb switch range start (nm) */
203 real rcoulomb
; /* Coulomb cutoff (nm) */
204 real epsilon_r
; /* relative dielectric constant */
205 real epsilon_rf
; /* relative dielectric constant of the RF */
206 int implicit_solvent
;/* No (=explicit water), or GBSA solvent models */
207 int gb_algorithm
; /* Algorithm to use for calculation Born radii */
208 int nstgbradii
; /* Frequency of updating Generalized Born radii */
209 real rgbradii
; /* Cutoff for GB radii calculation */
210 real gb_saltconc
; /* Salt concentration (M) for GBSA models */
211 real gb_epsilon_solvent
; /* dielectric coeff. of implicit solvent */
212 real gb_obc_alpha
; /* 1st scaling factor for Bashford-Case GB */
213 real gb_obc_beta
; /* 2nd scaling factor for Bashford-Case GB */
214 real gb_obc_gamma
; /* 3rd scaling factor for Bashford-Case GB */
215 real gb_dielectric_offset
; /* Dielectric offset for Still/HCT/OBC */
216 int sa_algorithm
; /* Algorithm for SA part of GBSA */
217 real sa_surface_tension
; /* Energy factor for SA part of GBSA */
218 int vdwtype
; /* Type of Van der Waals treatment */
219 real rvdw_switch
; /* Van der Waals switch range start (nm) */
220 real rvdw
; /* Van der Waals cutoff (nm) */
221 int eDispCorr
; /* Perform Long range dispersion corrections */
222 real tabext
; /* Extension of the table beyond the cut-off, *
223 * as well as the table length for 1-4 interac. */
224 real shake_tol
; /* tolerance for shake */
225 int efep
; /* free energy interpolation no/yes */
226 double init_lambda
; /* initial value for perturbation variable */
227 double delta_lambda
; /* change of lambda per time step (1/dt) */
228 int n_flambda
; /* The number of foreign lambda points */
229 double *flambda
; /* The foreign lambda values */
230 real sc_alpha
; /* free energy soft-core parameter */
231 int sc_power
; /* lambda power for soft-core interactions */
232 real sc_sigma
; /* free energy soft-core sigma when c6 or c12=0 */
233 real sc_sigma_min
; /* minimum FE sc sigma (default: =sg_sigma) */
234 int nstdhdl
; /* The frequency for writing to dhdl.xvg */
235 int dh_table_size
; /* The maximum table size for the dH table */
236 double dh_table_spacing
; /* The spacing for the dH table */
237 int eDisre
; /* Type of distance restraining */
238 real dr_fc
; /* force constant for ta_disre */
239 int eDisreWeighting
; /* type of weighting of pairs in one restraints */
240 bool bDisreMixed
; /* Use comb of time averaged and instan. viol's */
241 int nstdisreout
; /* frequency of writing pair distances to enx */
242 real dr_tau
; /* time constant for memory function in disres */
243 real orires_fc
; /* force constant for orientational restraints */
244 real orires_tau
; /* time constant for memory function in orires */
245 int nstorireout
; /* frequency of writing tr(SD) to enx */
246 real dihre_fc
; /* force constant for dihedral restraints */
247 real em_stepsize
; /* The stepsize for updating */
248 real em_tol
; /* The tolerance */
249 int niter
; /* Number of iterations for convergence of */
250 /* steepest descent in relax_shells */
251 real fc_stepsize
; /* Stepsize for directional minimization */
252 /* in relax_shells */
253 int nstcgsteep
; /* number of steps after which a steepest */
254 /* descents step is done while doing cg */
255 int nbfgscorr
; /* Number of corrections to the hessian to keep */
256 int eConstrAlg
; /* Type of constraint algorithm */
257 int nProjOrder
; /* Order of the LINCS Projection Algorithm */
258 real LincsWarnAngle
; /* If bond rotates more than %g degrees, warn */
259 int nLincsIter
; /* Number of iterations in the final Lincs step */
260 bool bShakeSOR
; /* Use successive overrelaxation for shake */
261 real bd_fric
; /* Friction coefficient for BD (amu/ps) */
262 int ld_seed
; /* Random seed for SD and BD */
263 int nwall
; /* The number of walls */
264 int wall_type
; /* The type of walls */
265 real wall_r_linpot
; /* The potentail is linear for r<=wall_r_linpot */
266 int wall_atomtype
[2];/* The atom type for walls */
267 real wall_density
[2]; /* Number density for walls */
268 real wall_ewald_zfac
; /* Scaling factor for the box for Ewald */
269 int ePull
; /* Type of pulling: no, umbrella or constraint */
270 t_pull
*pull
; /* The data for center of mass pulling */
271 real cos_accel
; /* Acceleration for viscosity calculation */
272 tensor deform
; /* Triclinic deformation velocities (nm/ps) */
273 int userint1
; /* User determined parameters */
281 t_grpopts opts
; /* Group options */
282 t_cosines ex
[DIM
]; /* Electric field stuff (spatial part) */
283 t_cosines et
[DIM
]; /* Electric field stuff (time part) */
284 bool bQMMM
; /* QM/MM calculation */
285 int QMconstraints
; /* constraints on QM bonds */
286 int QMMMscheme
; /* Scheme: ONIOM or normal */
287 real scalefactor
; /* factor for scaling the MM charges in QM calc.*/
290 #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)
292 #define DYNAMIC_BOX(ir) ((ir).epc!=epcNO || (ir).eI==eiTPI || DEFORM(ir))
294 #define PRESERVE_SHAPE(ir) ((ir).epc != epcNO && (ir).deform[XX][XX] == 0 && ((ir).epct == epctISOTROPIC || (ir).epct == epctSEMIISOTROPIC))
296 #define NEED_MUTOT(ir) (((ir).coulombtype==eelEWALD || EEL_PME((ir).coulombtype)) && ((ir).ewald_geometry==eewg3DC || (ir).epsilon_surface!=0))
298 #define IR_TWINRANGE(ir) ((ir).rlist > 0 && ((ir).rlistlong == 0 || (ir).rlistlong > (ir).rlist))
300 #define IR_ELEC_FIELD(ir) ((ir).ex[XX].n > 0 || (ir).ex[YY].n > 0 || (ir).ex[ZZ].n > 0)
302 #define IR_EXCL_FORCES(ir) (EEL_FULL((ir).coulombtype) || (EEL_RF((ir).coulombtype) && (ir).coulombtype != eelRF_NEC) || (ir).implicit_solvent != eisNO)
303 /* use pointer definitions of ir here, since that's what's usually used in the code */
304 #define IR_NVT_TROTTER(ir) ((((ir)->eI == eiVV) || ((ir)->eI == eiVVAK)) && ((ir)->etc == etcNOSEHOOVER))
306 #define IR_NPT_TROTTER(ir) ((((ir)->eI == eiVV) || ((ir)->eI == eiVVAK)) && ((ir)->epc == epcMTTK))