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