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
40 #include "../sysstuff.h"
48 int n; /* Number of terms */
49 real *a; /* Coeffients (V / nm ) */
50 real *phi; /* Phase angles */
54 real E0; /* Field strength (V/nm) */
55 real omega; /* Frequency (1/ps) */
56 real t0; /* Centre of the Gaussian pulse (ps) */
57 real sigma; /* Width of the Gaussian pulse (FWHM) (ps) */
60 #define EGP_EXCL (1<<0)
61 #define EGP_TABLE (1<<1)
64 int ngtc; /* # T-Coupl groups */
65 int nhchainlength; /* # of nose-hoover chains per group */
66 int ngacc; /* # Accelerate groups */
67 int ngfrz; /* # Freeze groups */
68 int ngener; /* # Ener groups */
69 real *nrdf; /* Nr of degrees of freedom in a group */
70 real *ref_t; /* Coupling temperature per group */
71 int *annealing; /* No/simple/periodic SA for each group */
72 int *anneal_npoints; /* Number of annealing time points per grp */
73 real **anneal_time; /* For ea. group: Time points */
74 real **anneal_temp; /* For ea. grp: Temperature at these times */
75 /* Final temp after all intervals is ref_t */
76 real *tau_t; /* Tau coupling time */
77 rvec *acc; /* Acceleration per group */
78 ivec *nFreeze; /* Freeze the group in each direction ? */
79 int *egp_flags; /* Exclusions/tables of energy group pairs */
82 int ngQM; /* nr of QM groups */
83 int *QMmethod; /* Level of theory in the QM calculation */
84 int *QMbasis; /* Basisset in the QM calculation */
85 int *QMcharge; /* Total charge in the QM region */
86 int *QMmult; /* Spin multiplicicty in the QM region */
87 gmx_bool *bSH; /* surface hopping (diabatic hop only) */
88 int *CASorbitals; /* number of orbiatls in the active space */
89 int *CASelectrons;/* number of electrons in the active space */
90 real *SAon; /* at which gap (A.U.) the SA is switched on */
92 int *SAsteps; /* in how many steps SA goes from 1-1 to 0.5-0.5*/
97 enum { epgrppbcNONE, epgrppbcREFAT, epgrppbcCOS };
100 int nat; /* Number of atoms in the pull group */
101 atom_id *ind; /* The global atoms numbers */
102 int nat_loc; /* Number of local pull atoms */
103 int nalloc_loc; /* Allocation size for ind_loc and weight_loc */
104 atom_id *ind_loc; /* Local pull indices */
105 int nweight; /* The number of weights (0 or nat) */
106 real *weight; /* Weights (use all 1 when weight==NULL) */
107 real *weight_loc; /* Weights for the local indices */
108 int epgrppbc; /* The type of pbc for this pull group, see enum above */
109 atom_id pbcatom; /* The reference atom for pbc (global number) */
110 rvec vec; /* The pull vector, direction or position */
111 rvec init; /* Initial reference displacement */
112 real rate; /* Rate of motion (nm/ps) */
113 real k; /* force constant */
114 real kB; /* force constant for state B */
115 real wscale; /* scaling factor for the weights: sum w m/sum w w m */
116 real invtm; /* inverse total mass of the group: 1/wscale sum w m */
117 dvec x; /* center of mass before update */
118 dvec xp; /* center of mass after update before constraining */
119 dvec dr; /* The distance from the reference group */
120 double f_scal; /* Scalar force for directional pulling */
121 dvec f; /* force due to the pulling/constraining */
125 int ngrp; /* number of groups */
126 int eGeom; /* pull geometry */
127 ivec dim; /* used to select components for constraint */
128 real cyl_r1; /* radius of cylinder for dynamic COM */
129 real cyl_r0; /* radius of cylinder including switch length */
130 real constr_tol; /* absolute tolerance for constraints in (nm) */
131 int nstxout; /* Output frequency for pull x */
132 int nstfout; /* Output frequency for pull f */
133 int ePBC; /* the boundary conditions */
134 int npbcdim; /* do pbc in dims 0 <= dim < npbcdim */
135 gmx_bool bRefAt; /* do we need reference atoms for a group COM ? */
136 int cosdim; /* dimension for cosine weighting, -1 if none */
137 gmx_bool bVirial; /* do we need to add the pull virial? */
138 t_pullgrp *grp; /* groups to pull/restrain/etc/ */
139 t_pullgrp *dyna; /* dynamic groups for use with local constraints */
140 rvec *rbuf; /* COM calculation buffer */
141 dvec *dbuf; /* COM calculation buffer */
142 double *dbuf_cyl; /* cylinder ref. groups COM calculation buffer */
144 FILE *out_x; /* output file for pull data */
145 FILE *out_f; /* output file for pull data */
149 /* Abstract types for enforced rotation only defined in pull_rotation.c */
150 typedef struct gmx_enfrot *gmx_enfrot_t;
151 typedef struct gmx_enfrotgrp *gmx_enfrotgrp_t;
154 int eType; /* Rotation type for this group */
155 int bMassW; /* Use mass-weighed positions? */
156 int nat; /* Number of atoms in the group */
157 atom_id *ind; /* The global atoms numbers */
158 rvec *x_ref; /* The reference positions */
159 rvec vec; /* The normalized rotation vector */
160 real rate; /* Rate of rotation (degree/ps) */
161 real k; /* Force constant (kJ/(mol nm^2) */
162 rvec pivot; /* Pivot point of rotation axis (nm) */
163 int eFittype; /* Type of fit to determine actual group angle */
164 real slab_dist; /* Slab distance (nm) */
165 real min_gaussian; /* Minimum value the gaussian must have so that
166 the force is actually evaluated */
167 real eps; /* Additive constant for radial motion2 and
168 flexible2 potentials (nm^2) */
169 gmx_enfrotgrp_t enfrotgrp; /* Stores non-inputrec rotation data per group */
173 int ngrp; /* Number of rotation groups */
174 int nstrout; /* Output frequency for main rotation outfile */
175 int nstsout; /* Output frequency for per-slab data */
176 t_rotgrp *grp; /* Groups to rotate */
177 gmx_enfrot_t enfrot; /* Stores non-inputrec enforced rotation data */
182 int eI; /* Integration method */
183 gmx_large_int_t nsteps; /* number of steps to be taken */
184 int simulation_part; /* Used in checkpointing to separate chunks */
185 gmx_large_int_t init_step; /* start at a stepcount >0 (used w. tpbconv) */
186 int nstcalcenergy; /* fequency of energy calc. and T/P coupl. upd. */
187 int ns_type; /* which ns method should we use? */
188 int nstlist; /* number of steps before pairlist is generated */
189 int ndelta; /* number of cells per rlong */
190 int nstcomm; /* number of steps after which center of mass */
191 /* motion is removed */
192 int comm_mode; /* Center of mass motion removal algorithm */
193 int nstcheckpoint; /* checkpointing frequency */
194 int nstlog; /* number of steps after which print to logfile */
195 int nstxout; /* number of steps after which X is output */
196 int nstvout; /* id. for V */
197 int nstfout; /* id. for F */
198 int nstenergy; /* number of steps after which energies printed */
199 int nstxtcout; /* id. for compressed trj (.xtc) */
200 double init_t; /* initial time (ps) */
201 double delta_t; /* time step (ps) */
202 real xtcprec; /* precision of xtc file */
203 int nkx,nky,nkz; /* number of k vectors in each spatial dimension*/
204 /* for fourier methods for long range electrost.*/
205 int pme_order; /* interpolation order for PME */
206 real ewald_rtol; /* Real space tolerance for Ewald, determines */
207 /* the real/reciprocal space relative weight */
208 int ewald_geometry; /* normal/3d ewald, or pseudo-2d LR corrections */
209 real epsilon_surface; /* Epsilon for PME dipole correction */
210 gmx_bool bOptFFT; /* optimize the fft plan at start */
211 int ePBC; /* Type of periodic boundary conditions */
212 int bPeriodicMols; /* Periodic molecules */
213 gmx_bool bContinuation; /* Continuation run: starting state is correct */
214 int etc; /* temperature coupling */
215 int nsttcouple; /* interval in steps for temperature coupling */
216 int epc; /* pressure coupling */
217 int epct; /* pressure coupling type */
218 int nstpcouple; /* interval in steps for pressure coupling */
219 real tau_p; /* pressure coupling time (ps) */
220 tensor ref_p; /* reference pressure (kJ/(mol nm^3)) */
221 tensor compress; /* compressability ((mol nm^3)/kJ) */
222 int refcoord_scaling;/* How to scale absolute reference coordinates */
223 rvec posres_com; /* The COM of the posres atoms */
224 rvec posres_comB; /* The B-state COM of the posres atoms */
225 int andersen_seed; /* Random seed for Andersen thermostat. */
226 real rlist; /* short range pairlist cut-off (nm) */
227 real rlistlong; /* long range pairlist cut-off (nm) */
228 real rtpi; /* Radius for test particle insertion */
229 int coulombtype; /* Type of electrostatics treatment */
230 real rcoulomb_switch; /* Coulomb switch range start (nm) */
231 real rcoulomb; /* Coulomb cutoff (nm) */
232 real epsilon_r; /* relative dielectric constant */
233 real epsilon_rf; /* relative dielectric constant of the RF */
234 int implicit_solvent;/* No (=explicit water), or GBSA solvent models */
235 int gb_algorithm; /* Algorithm to use for calculation Born radii */
236 int nstgbradii; /* Frequency of updating Generalized Born radii */
237 real rgbradii; /* Cutoff for GB radii calculation */
238 real gb_saltconc; /* Salt concentration (M) for GBSA models */
239 real gb_epsilon_solvent; /* dielectric coeff. of implicit solvent */
240 real gb_obc_alpha; /* 1st scaling factor for Bashford-Case GB */
241 real gb_obc_beta; /* 2nd scaling factor for Bashford-Case GB */
242 real gb_obc_gamma; /* 3rd scaling factor for Bashford-Case GB */
243 real gb_dielectric_offset; /* Dielectric offset for Still/HCT/OBC */
244 int sa_algorithm; /* Algorithm for SA part of GBSA */
245 real sa_surface_tension; /* Energy factor for SA part of GBSA */
246 int vdwtype; /* Type of Van der Waals treatment */
247 real rvdw_switch; /* Van der Waals switch range start (nm) */
248 real rvdw; /* Van der Waals cutoff (nm) */
249 int eDispCorr; /* Perform Long range dispersion corrections */
250 real tabext; /* Extension of the table beyond the cut-off, *
251 * as well as the table length for 1-4 interac. */
252 real shake_tol; /* tolerance for shake */
253 int efep; /* free energy interpolation no/yes */
254 double init_lambda; /* initial value for perturbation variable */
255 double delta_lambda; /* change of lambda per time step (1/dt) */
256 int n_flambda; /* The number of foreign lambda points */
257 double *flambda; /* The foreign lambda values */
258 real sc_alpha; /* free energy soft-core parameter */
259 int sc_power; /* lambda power for soft-core interactions */
260 real sc_sigma; /* free energy soft-core sigma when c6 or c12=0 */
261 real sc_sigma_min; /* minimum FE sc sigma (default: =sg_sigma) */
262 int nstdhdl; /* The frequency for writing to dhdl.xvg */
263 int separate_dhdl_file; /* whether to write a separate dhdl.xvg file
264 note: NOT a gmx_bool, but an enum */
265 int dhdl_derivatives;/* whether to calculate+write dhdl derivatives
266 note: NOT a gmx_bool, but an enum */
267 int dh_hist_size; /* The maximum size for the dH histogram */
268 double dh_hist_spacing; /* The spacing for the dH histogram */
269 int eDisre; /* Type of distance restraining */
270 real dr_fc; /* force constant for ta_disre */
271 int eDisreWeighting; /* type of weighting of pairs in one restraints */
272 gmx_bool bDisreMixed; /* Use comb of time averaged and instan. viol's */
273 int nstdisreout; /* frequency of writing pair distances to enx */
274 real dr_tau; /* time constant for memory function in disres */
275 real orires_fc; /* force constant for orientational restraints */
276 real orires_tau; /* time constant for memory function in orires */
277 int nstorireout; /* frequency of writing tr(SD) to enx */
278 real dihre_fc; /* force constant for dihedral restraints */
279 real em_stepsize; /* The stepsize for updating */
280 real em_tol; /* The tolerance */
281 int niter; /* Number of iterations for convergence of */
282 /* steepest descent in relax_shells */
283 real fc_stepsize; /* Stepsize for directional minimization */
284 /* in relax_shells */
285 int nstcgsteep; /* number of steps after which a steepest */
286 /* descents step is done while doing cg */
287 int nbfgscorr; /* Number of corrections to the hessian to keep */
288 int eConstrAlg; /* Type of constraint algorithm */
289 int nProjOrder; /* Order of the LINCS Projection Algorithm */
290 real LincsWarnAngle; /* If bond rotates more than %g degrees, warn */
291 int nLincsIter; /* Number of iterations in the final Lincs step */
292 gmx_bool bShakeSOR; /* Use successive overrelaxation for shake */
293 real bd_fric; /* Friction coefficient for BD (amu/ps) */
294 int ld_seed; /* Random seed for SD and BD */
295 int nwall; /* The number of walls */
296 int wall_type; /* The type of walls */
297 real wall_r_linpot; /* The potentail is linear for r<=wall_r_linpot */
298 int wall_atomtype[2];/* The atom type for walls */
299 real wall_density[2]; /* Number density for walls */
300 real wall_ewald_zfac; /* Scaling factor for the box for Ewald */
301 int ePull; /* Type of pulling: no, umbrella or constraint */
302 t_pull *pull; /* The data for center of mass pulling */
303 gmx_bool bRot; /* Calculate enforced rotation potential(s)? */
304 t_rot *rot; /* The data for enforced rotation potentials */
305 real cos_accel; /* Acceleration for viscosity calculation */
306 tensor deform; /* Triclinic deformation velocities (nm/ps) */
307 int userint1; /* User determined parameters */
315 t_grpopts opts; /* Group options */
316 t_cosines ex[DIM]; /* Electric field stuff (spatial part) */
317 t_cosines et[DIM]; /* Electric field stuff (time part) */
318 gmx_bool bQMMM; /* QM/MM calculation */
319 int QMconstraints; /* constraints on QM bonds */
320 int QMMMscheme; /* Scheme: ONIOM or normal */
321 real scalefactor; /* factor for scaling the MM charges in QM calc.*/
324 #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)
326 #define DYNAMIC_BOX(ir) ((ir).epc!=epcNO || (ir).eI==eiTPI || DEFORM(ir))
328 #define PRESERVE_SHAPE(ir) ((ir).epc != epcNO && (ir).deform[XX][XX] == 0 && ((ir).epct == epctISOTROPIC || (ir).epct == epctSEMIISOTROPIC))
330 #define NEED_MUTOT(ir) (((ir).coulombtype==eelEWALD || EEL_PME((ir).coulombtype)) && ((ir).ewald_geometry==eewg3DC || (ir).epsilon_surface!=0))
332 #define IR_TWINRANGE(ir) ((ir).rlist > 0 && ((ir).rlistlong == 0 || (ir).rlistlong > (ir).rlist))
334 #define IR_ELEC_FIELD(ir) ((ir).ex[XX].n > 0 || (ir).ex[YY].n > 0 || (ir).ex[ZZ].n > 0)
336 #define IR_EXCL_FORCES(ir) (EEL_FULL((ir).coulombtype) || (EEL_RF((ir).coulombtype) && (ir).coulombtype != eelRF_NEC) || (ir).implicit_solvent != eisNO)
337 /* use pointer definitions of ir here, since that's what's usually used in the code */
338 #define IR_NVT_TROTTER(ir) ((((ir)->eI == eiVV) || ((ir)->eI == eiVVAK)) && ((ir)->etc == etcNOSEHOOVER))
340 #define IR_NPT_TROTTER(ir) ((((ir)->eI == eiVV) || ((ir)->eI == eiVVAK)) && ((ir)->epc == epcMTTK))