4 * This source code is part of
8 * GROningen MAchine for Chemical Simulations
12 * Copyright (c) 1991-1999
13 * BIOSON Research Institute, Dept. of Biophysical Chemistry
14 * University of Groningen, The Netherlands
17 * GROMACS: A message-passing parallel molecular dynamics implementation
18 * H.J.C. Berendsen, D. van der Spoel and R. van Drunen
19 * Comp. Phys. Comm. 91, 43-56 (1995)
21 * Also check out our WWW page:
22 * http://md.chem.rug.nl/~gmx
27 * Good ROcking Metal Altar for Chronical Sinners
33 static char *SRCID_gstat_h
= "$Id$";
44 /***********************************************
46 * A U T O C O R R E L A T I O N
48 ***********************************************/
50 extern real
LegendreP(real x
,unsigned long m
);
52 #define eacNormal (1<<0)
54 #define eacVector (1<<2)
55 #define eacRcross (1<<3 | eacVector)
56 #define eacP0 (1<<4 | eacVector)
57 #define eacP1 (1<<5 | eacVector)
58 #define eacP2 (1<<6 | eacVector)
59 #define eacP3 (1<<7 | eacVector)
60 #define eacP4 (1<<8 | eacVector)
62 extern t_pargs
*add_acf_pargs(int *npargs
,t_pargs
*pa
);
63 /* Add options for autocorr to the current set of options.
64 * *npargs must be initialised to the number of elements in pa,
65 * it will be incremented appropriately.
68 extern void do_autocorr(char *fn
,char *title
,int nframes
,int nitem
,real
**c1
,
69 real dt
,unsigned long mode
,bool bAver
);
70 /* Calls low_do_autocorr (see below). After calling add_acf_pargs */
72 extern void low_do_autocorr(char *fn
,char *title
,
73 int nframes
,int nitem
,int nout
,real
**c1
,
74 real dt
,unsigned long mode
,int nrestart
,
75 bool bAver
,bool bFour
,bool bNormalize
,
76 bool bVerbose
,real tbeginfit
,real tendfit
,
79 * do_autocorr calculates autocorrelation functions for many things.
80 * It takes a 2 d array containing nitem arrays of length nframes
81 * for each item the ACF is calculated.
83 * A number of "modes" exist for computation of the ACF
85 * if (mode == eacNormal) {
86 * C(t) = < X (tau) * X (tau+t) >
88 * else if (mode == eacCos) {
89 * C(t) = < cos (X(tau) - X(tau+t)) >
91 * else if (mode == eacVector) {
92 * C(t) = < X(tau) * X(tau+t)
94 * else if (mode == eacP1) {
95 * C(t) = < cos (X(tau) * X(tau+t) >
97 * else if (mode == eacP2) {
98 * C(t) = 1/2 * < 3 cos (X(tau) * X(tau+t) - 1 >
100 * else if (mode == eacRcross) {
101 * C(t) = < ( X(tau) * X(tau+t) )^2 >
104 * For modes eacVector, eacP1, eacP2 and eacRcross the input should be
105 * 3 x nframes long, where each triplet is taken as a 3D vector
107 * For mode eacCos inputdata must be in radians, not degrees!
109 * Other parameters are:
111 * fn is output filename (.xvg) where the correlation function(s) are printed
112 * title is the title in the output file
113 * nframes is the number of frames in the time series
114 * nitem is the number of items
115 * c1 is an array of dimension [ 0 .. nitem-1 ] [ 0 .. nframes-1 ]
116 * on output, this array is filled with the correlation function
118 * nrestart is the number of steps between restarts for direct ACFs
119 * (i.e. without FFT) When set to 1 all points are used as
120 * time origin for averaging
121 * dt is the time between frames
122 * bAver If set, all ndih C(t) functions are averaged into a single
124 * bFour If set, will use fast fourier transform (FFT) for evaluating
126 * bNormalize If set, all ACFs will be normalized to start at 0
130 char *name
; /* Description of the J coupling constant */
131 real A
,B
,C
; /* Karplus coefficients */
132 real offset
; /* Offset for dihedral angle in histogram (e.g. -M_PI/3) */
133 real Jc
; /* Resulting Jcoupling */
136 extern void calc_distribution_props(int nh
,int histo
[],
137 real start
,int nkkk
, t_karplus kkk
[],
139 /* This routine takes a dihedral distribution and calculates
140 * coupling constants and dihedral order parameters of it.
142 * nh is the number of points
143 * histo is the array of datapoints which is assumed to span
145 * start is the starting angle of the histogram, this can be either 0
147 * nkkk is the number of karplus sets (multiple coupling constants may be
148 * derived from a single angle)
149 * kkk are the constants for calculating J coupling constants using a
150 * Karplus equation according to
153 * J = A cos theta + B cos theta + C
155 * where theta is phi - offset (phi is the angle in the histogram)
156 * offset is subtracted from phi before substitution in the Karplus
158 * S2 is the resulting dihedral order parameter
163 /***********************************************
165 * F I T R O U T I N E S
167 ***********************************************/
168 extern void do_expfit(int ndata
,real c1
[],real dt
,
169 real begintimefit
,real endtimefit
);
171 extern void expfit(int n
, real x
[], real y
[], real Dy
[],
174 /* This procedure fits y=exp(a+bx) for n (x,y) pairs to determine a and b.
175 * The uncertainties in the y values must be in the vector Dy.
176 * The standard deviations of a and b, sa and sb, are also calculated.
178 * Routine from Computers in physics, 7(3) (1993), p. 280-285.
181 extern void ana_dih_trans(char *fn_trans
,char *fn_histo
,
182 real
**dih
,int nframes
,int nangles
,
183 char *grpname
,real t0
,real dt
,bool bRb
);
185 * Analyse dihedral transitions, by counting transitions per dihedral
186 * and per frame. The total number of transitions is printed to
187 * stderr, as well as the average time between transitions.
189 * Dihedrals are supposed to be in either of three minima,
190 * (trans, gauche+, gauche-)
192 * fn_trans output file name for #transitions per timeframe
193 * fn_histo output file name for transition time histogram
194 * dih the actual dihedral angles
195 * nframes number of times frames
196 * nangles number of angles
197 * grpname a string for the header of plots
198 * t0,dt starting time resp. time between time frames
199 * bRb determines whether the polymer convention is used
203 extern void read_ang_dih(char *trj_fn
,char *tpb_fn
,
204 bool bAngles
,bool bSaveAll
,bool bRb
,
205 int maxangstat
,int angstat
[],
206 int *nframes
,real
**time
,
207 int isize
,atom_id index
[],
212 * Read a trajectory and calculate angles and dihedrals.
214 * trj_fn file name of trajectory
215 * tpb_fn file name of tpb file
216 * bAngles do we have to read angles or dihedrals
217 * bSaveAll do we have to store all in the dih array
218 * bRb do we have Ryckaert-Bellemans dihedrals (trans = 0)
219 * maxangstat number of entries in distribution array
220 * angstat angle distribution
221 * *nframes number of frames read
222 * time simulation time at each time frame
223 * isize number of entries in the index, when angles 3*number of angles
224 * else 4*number of angles
225 * index atom numbers that define the angles or dihedrals
226 * (i,j,k) resp (i,j,k,l)
227 * trans_frac number of dihedrals in trans
228 * aver_angle average angle at each time frame
229 * dih all angles at each time frame
232 extern void make_histo(FILE *log
,
233 int ndata
,real data
[],int npoints
,int histo
[],
234 real minx
,real maxx
);
236 * Make a histogram from data. The min and max of the data array can
237 * be determined (if minx == 0 and maxx == 0)
238 * and the index in the histogram is computed from
239 * ind = npoints/(max(data) - min(data))
241 * log write error output to this file
242 * ndata number of points in data
244 * npoints number of points in histogram
245 * histo histogram array. This is NOT set to zero, to allow you
246 * to add multiple histograms
247 * minx start of the histogram
248 * maxx end of the histogram
249 * if both are 0, these values are computed by the routine itself
252 extern void normalize_histo(int npoints
,int histo
[],real dx
,real normhisto
[]);
254 * Normalize a histogram so that the integral over the histo is 1
256 * npoints number of points in the histo array
257 * histo input histogram
258 * dx distance between points on the X-axis
259 * normhisto normalized output histogram
262 /* Use Levenberg-Marquardt method to fit to a one parameter exponential */
263 /* Or: "There is no KILL like OVERKILL", Dr. Ir. D. van der Spoel */
264 extern real
do_lmfit(int ndata
,real c1
[],real sig
[],real dt
,
265 real begintimefit
,real endtimefit
,
266 bool bVerbose
,int nfitparm
,
267 real fit
[],real fitparms
[],char *fix
);
268 /* Returns integral. if fit != NULL, than the fitted function is copied
269 * into the original data array
272 extern real
print_and_integrate(FILE *fp
,int n
,real dt
,real c
[]);
273 /* Integrate the data in c[] from 0 to n using trapezium rule.
274 * If fp != NULL output is written to it
277 extern int get_acfnout(void);
278 /* Return the output length for the correlation function
279 * Works only AFTER do_auto_corr has been called!
283 /* Least squares method of computing statistics */
285 double yy
,yx
,xx
,sx
,sy
;
289 extern void init_lsq(t_lsq
*lsq
);
290 extern void done_lsq(t_lsq
*lsq
);
291 extern void add_lsq(t_lsq
*lsq
,real x
,real y
);
292 extern void get_lsq_ab(t_lsq
*lsq
,real
*a
,real
*b
);
293 extern int npoints_lsq(t_lsq
*lsq
);
294 extern real
aver_lsq(t_lsq
*lsq
);
295 extern real
sigma_lsq(t_lsq
*lsq
);
296 extern real
error_lsq(t_lsq
*lsq
);