include/gstat.h

   1 /*
   2  * $Id$
   3  *
   4  *                This source code is part of
   5  *
   6  *                 G   R   O   M   A   C   S
   7  *
   8  *          GROningen MAchine for Chemical Simulations
   9  *
  10  *                        VERSION 3.0
  11  *
  12  * Copyright (c) 1991-2001
  13  * BIOSON Research Institute, Dept. of Biophysical Chemistry
  14  * University of Groningen, The Netherlands
  15  *
  16  * This program is free software; you can redistribute it and/or
  17  * modify it under the terms of the GNU General Public License
  18  * as published by the Free Software Foundation; either version 2
  19  * of the License, or (at your option) any later version.
  20  *
  21  * If you want to redistribute modifications, please consider that
  22  * scientific software is very special. Version control is crucial -
  23  * bugs must be traceable. We will be happy to consider code for
  24  * inclusion in the official distribution, but derived work must not
  25  * be called official GROMACS. Details are found in the README & COPYING
  26  * files - if they are missing, get the official version at www.gromacs.org.
  27  *
  28  * To help us fund GROMACS development, we humbly ask that you cite
  29  * the papers on the package - you can find them in the top README file.
  30  *
  31  * Do check out http://www.gromacs.org , or mail us at gromacs@gromacs.org .
  32  *
  33  * And Hey:
  34  * GRoups of Organic Molecules in ACtion for Science
  35  */
  36
  37 #ifndef _gstat_h
  38 #define _gstat_h
  39
  40 static char *SRCID_gstat_h = "$Id$";
  41 #ifdef HAVE_CONFIG_H
  42 #include <config.h>
  43 #endif
  44
  45 static char *SRCID_gstat_h = "$Id$";
  46
  47 #include "typedefs.h"
  48 #include "statutil.h"
  49 #include "mshift.h"
  50 #include "rmpbc.h"
  51
  52 #ifdef CPLUSPLUS
  53 extern "C" {
  54 #endif
  55
  56 /***********************************************
  57  *
  58  *     A U T O C O R R E L A T I O N
  59  *
  60  ***********************************************/
  61
  62 extern real LegendreP(real x,unsigned long m);
  63
  64 #define eacNormal (1<<0)
  65 #define eacCos    (1<<1)
  66 #define eacVector (1<<2)
  67 #define eacRcross (1<<3  | eacVector)
  68 #define eacP0     (1<<4  | eacVector)
  69 #define eacP1     (1<<5  | eacVector)
  70 #define eacP2     (1<<6  | eacVector)
  71 #define eacP3     (1<<7  | eacVector)
  72 #define eacP4     (1<<8  | eacVector)
  73
  74 enum {
  75   effnNONE, effnEXP1, effnEXP2, effnEXP3, effnVAC, effnERREST, effnNR
  76 };
  77
  78 extern int  nfp_ffn[effnNR];
  79
  80 extern char *s_ffn[effnNR+2];
  81
  82 extern char *longs_ffn[effnNR];
  83
  84 int sffn2effn(char **sffn);
  85 /* Returns the ffn enum corresponding to the selected enum option in sffn */
  86
  87 extern t_pargs *add_acf_pargs(int *npargs,t_pargs *pa);
  88 /* Add options for autocorr to the current set of options.
  89  * *npargs must be initialised to the number of elements in pa,
  90  * it will be incremented appropriately.
  91  */
  92
  93 extern void do_autocorr(char *fn,char *title,int nframes,int nitem,real **c1,
  94                         real dt,unsigned long mode,bool bAver);
  95 /* Calls low_do_autocorr (see below). After calling add_acf_pargs */
  96
  97 extern void low_do_autocorr(char *fn,char *title,
  98                             int  nframes,int nitem,int nout,real **c1,
  99                             real dt,unsigned long mode,int nrestart,
 100                             bool bAver,bool bFour,bool bNormalize,
 101                             bool bVerbose,real tbeginfit,real tendfit,
 102                             int nfitparm,int nskip);
 103 /*
 104  * do_autocorr calculates autocorrelation functions for many things.
 105  * It takes a 2 d array containing nitem arrays of length nframes
 106  * for each item the ACF is calculated.
 107  *
 108  * A number of "modes" exist for computation of the ACF
 109  *
 110  * if (mode == eacNormal) {
 111  *   C(t) = < X (tau) * X (tau+t) >
 112  * }
 113  * else if (mode == eacCos) {
 114  *   C(t) = < cos (X(tau) - X(tau+t)) >
 115  * }
 116  * else if (mode == eacVector) {
 117  *   C(t) = < X(tau) * X(tau+t)
 118  * }
 119  * else if (mode == eacP1) {
 120  *   C(t) = < cos (X(tau) * X(tau+t) >
 121  * }
 122  * else if (mode == eacP2) {
 123  *   C(t) = 1/2 * < 3 cos (X(tau) * X(tau+t) - 1 >
 124  * }
 125  * else if (mode == eacRcross) {
 126  *   C(t) = < ( X(tau) * X(tau+t) )^2 >
 127  * }
 128  *
 129  * For modes eacVector, eacP1, eacP2 and eacRcross the input should be
 130  * 3 x nframes long, where each triplet is taken as a 3D vector
 131  *
 132  * For mode eacCos inputdata must be in radians, not degrees!
 133  *
 134  * Other parameters are:
 135  *
 136  * fn is output filename (.xvg) where the correlation function(s) are printed
 137  * title is the title in the output file
 138  * nframes is the number of frames in the time series
 139  * nitem is the number of items
 140  * c1           is an array of dimension [ 0 .. nitem-1 ] [ 0 .. nframes-1 ]
 141  *              on output, this array is filled with the correlation function
 142  *              to reduce storage
 143  * nrestart     is the number of steps between restarts for direct ACFs
 144  *              (i.e. without FFT) When set to 1 all points are used as
 145  *              time origin for averaging
 146  * dt           is the time between frames
 147  * bAver        If set, all ndih C(t) functions are averaged into a single
 148  *              C(t)
 149  * bFour        If set, will use fast fourier transform (FFT) for evaluating
 150  *              the ACF
 151  * bNormalize   If set, all ACFs will be normalized to start at 0
 152  * nskip        Determines whether steps a re skipped in the output
 153  */
 154
 155 typedef struct {
 156   char *name;   /* Description of the J coupling constant */
 157   real A,B,C;   /* Karplus coefficients */
 158   real offset;  /* Offset for dihedral angle in histogram (e.g. -M_PI/3) */
 159   real Jc;      /* Resulting Jcoupling */
 160 } t_karplus;
 161
 162 extern void calc_distribution_props(int nh,int histo[],
 163                                     real start,int  nkkk, t_karplus kkk[],
 164                                     real *S2);
 165 /* This routine takes a dihedral distribution and calculates
 166  * coupling constants and dihedral order parameters of it.
 167  *
 168  * nh      is the number of points
 169  * histo   is the array of datapoints which is assumed to span
 170  *         2 M_PI radians
 171  * start   is the starting angle of the histogram, this can be either 0
 172  *         or -M_PI
 173  * nkkk    is the number of karplus sets (multiple coupling constants may be
 174  *         derived from a single angle)
 175  * kkk     are the constants for calculating J coupling constants using a
 176  *         Karplus equation according to
 177  *
 178  *                  2
 179  *         J = A cos theta + B cos theta + C
 180  *
 181  *         where theta is phi - offset (phi is the angle in the histogram)
 182  * offset  is subtracted from phi before substitution in the Karplus
 183  *         equation
 184  * S2      is the resulting dihedral order parameter
 185  *
 186  */
 187
 188
 189 /***********************************************
 190  *
 191  *     F I T   R O U T I N E S
 192  *
 193  ***********************************************/
 194 extern void do_expfit(int ndata,real c1[],real dt,
 195                       real begintimefit,real endtimefit);
 196
 197 extern void expfit(int n, real x[], real y[], real Dy[],
 198                    real *a, real *sa,
 199                    real *b, real *sb);
 200 /* This procedure fits y=exp(a+bx) for n (x,y) pairs to determine a and b.
 201  * The uncertainties in the y values must be in the vector Dy.
 202  * The standard deviations of a and b, sa and sb, are also calculated.
 203  *
 204  * Routine from Computers in physics, 7(3) (1993), p. 280-285.
 205  */
 206
 207 extern void ana_dih_trans(char *fn_trans,char *fn_histo,
 208                           real **dih,int nframes,int nangles,
 209                           char *grpname,real t0,real dt,bool bRb);
 210 /*
 211  * Analyse dihedral transitions, by counting transitions per dihedral
 212  * and per frame. The total number of transitions is printed to
 213  * stderr, as well as the average time between transitions.
 214  *
 215  * Dihedrals are supposed to be in either of three minima,
 216  * (trans, gauche+, gauche-)
 217  *
 218  * fn_trans  output file name for #transitions per timeframe
 219  * fn_histo  output file name for transition time histogram
 220  * dih       the actual dihedral angles
 221  * nframes   number of times frames
 222  * nangles   number of angles
 223  * grpname   a string for the header of plots
 224  * t0,dt     starting time resp. time between time frames
 225  * bRb       determines whether the polymer convention is used
 226  *           (trans = 0)
 227  */
 228
 229 extern void read_ang_dih(char *trj_fn,char *tpb_fn,
 230                          bool bAngles,bool bSaveAll,bool bRb,
 231                          int maxangstat,int angstat[],
 232                          int *nframes,real **time,
 233                          int isize,atom_id index[],
 234                          real **trans_frac,
 235                          real **aver_angle,
 236                          real *dih[]);
 237 /*
 238  * Read a trajectory and calculate angles and dihedrals.
 239  *
 240  * trj_fn      file name of trajectory
 241  * tpb_fn      file name of tpb file
 242  * bAngles     do we have to read angles or dihedrals
 243  * bSaveAll    do we have to store all in the dih array
 244  * bRb         do we have Ryckaert-Bellemans dihedrals (trans = 0)
 245  * maxangstat  number of entries in distribution array
 246  * angstat     angle distribution
 247  * *nframes    number of frames read
 248  * time        simulation time at each time frame
 249  * isize       number of entries in the index, when angles 3*number of angles
 250  *             else 4*number of angles
 251  * index       atom numbers that define the angles or dihedrals
 252  *             (i,j,k) resp (i,j,k,l)
 253  * trans_frac  number of dihedrals in trans
 254  * aver_angle  average angle at each time frame
 255  * dih         all angles at each time frame
 256  */
 257
 258 extern void make_histo(FILE *log,
 259                        int ndata,real data[],int npoints,int histo[],
 260                        real minx,real maxx);
 261 /*
 262  * Make a histogram from data. The min and max of the data array can
 263  * be determined (if minx == 0 and maxx == 0)
 264  * and the index in the histogram is computed from
 265  * ind = npoints/(max(data) - min(data))
 266  *
 267  * log       write error output to this file
 268  * ndata     number of points in data
 269  * data      data points
 270  * npoints   number of points in histogram
 271  * histo     histogram array. This is NOT set to zero, to allow you
 272  *           to add multiple histograms
 273  * minx      start of the histogram
 274  * maxx      end of the histogram
 275  *           if both are 0, these values are computed by the routine itself
 276  */
 277
 278 extern void normalize_histo(int npoints,int histo[],real dx,real normhisto[]);
 279 /*
 280  * Normalize a histogram so that the integral over the histo is 1
 281  *
 282  * npoints    number of points in the histo array
 283  * histo      input histogram
 284  * dx         distance between points on the X-axis
 285  * normhisto  normalized output histogram
 286  */
 287
 288 extern real fit_function(int eFitFn,real *parm,real x);
 289 /* Returns the value of fit function eFitFn at x */
 290
 291 /* Use Levenberg-Marquardt method to fit to a nfitparm parameter exponential */
 292 /* or to a transverse current autocorrelation function */
 293 /* Or: "There is no KILL like OVERKILL", Dr. Ir. D. van der Spoel */
 294 extern real do_lmfit(int ndata,real c1[],real sig[],real dt,real *x,
 295                      real begintimefit,real endtimefit,bool bVerbose,
 296                      int eFitFn,real fitparms[],int fix);
 297 /* Returns integral.
 298  * If x == NULL, the timestep dt will be used to create a time axis.
 299  * fix fixes fit parameter i at it's starting value, when the i'th bit
 300  * of fix is set.
 301  */
 302
 303 extern real print_and_integrate(FILE *fp,int n,real dt,
 304                                 real c[],real *fit,int nskip);
 305 /* Integrate the data in c[] from 0 to n using trapezium rule.
 306  * If fp != NULL output is written to it
 307  * nskip determines whether all elements are written to the output file
 308  * (written when i % nskip == 0)
 309  * If fit != NULL the fit is also written.
 310  */
 311
 312 extern int get_acfnout(void);
 313 /* Return the output length for the correlation function
 314  * Works only AFTER do_auto_corr has been called!
 315  */
 316
 317
 318 /* Least squares method of computing statistics */
 319 typedef struct {
 320   double yy,yx,xx,sx,sy;
 321   int    np;
 322 } t_lsq;
 323
 324 extern void init_lsq(t_lsq *lsq);
 325 extern void done_lsq(t_lsq *lsq);
 326 extern void add_lsq(t_lsq *lsq,real x,real y);
 327 extern void get_lsq_ab(t_lsq *lsq,real *a,real *b);
 328 extern int  npoints_lsq(t_lsq *lsq);
 329 extern real aver_lsq(t_lsq *lsq);
 330 extern real sigma_lsq(t_lsq *lsq);
 331 extern real error_lsq(t_lsq *lsq);
 332
 333 #ifdef CPLUSPLUS
 334              }
 335 #endif
 336
 337 #endif
 338
 339