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.1
  11  * Copyright (c) 1991-2001, University of Groningen, The Netherlands
  12  * This program is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU General Public License
  14  * as published by the Free Software Foundation; either version 2
  15  * of the License, or (at your option) any later version.
  16  *
  17  * If you want to redistribute modifications, please consider that
  18  * scientific software is very special. Version control is crucial -
  19  * bugs must be traceable. We will be happy to consider code for
  20  * inclusion in the official distribution, but derived work must not
  21  * be called official GROMACS. Details are found in the README & COPYING
  22  * files - if they are missing, get the official version at www.gromacs.org.
  23  *
  24  * To help us fund GROMACS development, we humbly ask that you cite
  25  * the papers on the package - you can find them in the top README file.
  26  *
  27  * For more info, check our website at http://www.gromacs.org
  28  *
  29  * And Hey:
  30  * Getting the Right Output Means no Artefacts in Calculating Stuff
  31  */
  32
  33 #ifndef _gstat_h
  34 #define _gstat_h
  35
  36
  37 #ifdef HAVE_CONFIG_H
  38 #include <config.h>
  39 #endif
  40
  41 #include "typedefs.h"
  42 #include "statutil.h"
  43 #include "mshift.h"
  44 #include "rmpbc.h"
  45
  46 #ifdef CPLUSPLUS
  47 extern "C" {
  48 #endif
  49
  50 /***********************************************
  51  *
  52  *     A U T O C O R R E L A T I O N
  53  *
  54  ***********************************************/
  55
  56 extern real LegendreP(real x,unsigned long m);
  57
  58 #define eacNormal (1<<0)
  59 #define eacCos    (1<<1)
  60 #define eacVector (1<<2)
  61 #define eacRcross (1<<3  | eacVector)
  62 #define eacP0     (1<<4  | eacVector)
  63 #define eacP1     (1<<5  | eacVector)
  64 #define eacP2     (1<<6  | eacVector)
  65 #define eacP3     (1<<7  | eacVector)
  66 #define eacP4     (1<<8  | eacVector)
  67
  68 enum {
  69   effnNONE, effnEXP1, effnEXP2, effnEXP3, effnVAC,
  70   effnEXP5, effnEXP7, effnERREST, effnNR
  71 };
  72
  73 /* must correspond with 'leg' g_chi.c:727 */
  74 enum { edPhi=0, edPsi, edOmega, edChi1, edChi2, edChi3, edChi4, edChi5, edChi6, edMax };
  75
  76 enum { edPrintST=0,edPrintRO } ;
  77
  78 #define NHISTO 360
  79 #define NONCHI 3
  80 #define MAXCHI edMax-NONCHI
  81 #define NROT 4  /* number of rotamers: 1=g(-), 2=t, 3=g(+), 0=other */
  82
  83 typedef struct {
  84   int minO,minC,H,N,C,O,Cn[MAXCHI+3];
  85 } t_dihatms; /* Cn[0]=N, Cn[1]=Ca, Cn[2]=Cb etc. */
  86
  87 typedef struct {
  88   char name[12];
  89   int  resnr;
  90   int  index;       /* Index for amino acids (histograms) */
  91   int  j0[edMax];   /* Index in dih array (phi angle is first...) */
  92   t_dihatms  atm;
  93   int  b[edMax];
  94   int  ntr[edMax];
  95   real S2[edMax];
  96   real rot_occ[edMax][NROT];
  97
  98 } t_dlist;
  99
 100 extern int  nfp_ffn[effnNR];
 101
 102 extern char *s_ffn[effnNR+2];
 103
 104 extern char *longs_ffn[effnNR];
 105
 106 int sffn2effn(char **sffn);
 107 /* Returns the ffn enum corresponding to the selected enum option in sffn */
 108
 109 extern t_pargs *add_acf_pargs(int *npargs,t_pargs *pa);
 110 /* Add options for autocorr to the current set of options.
 111  * *npargs must be initialised to the number of elements in pa,
 112  * it will be incremented appropriately.
 113  */
 114
 115 extern void do_autocorr(char *fn,char *title,int nframes,int nitem,real **c1,
 116                         real dt,unsigned long mode,bool bAver);
 117 /* Calls low_do_autocorr (see below). After calling add_acf_pargs */
 118
 119 extern void low_do_autocorr(char *fn,char *title,
 120                             int  nframes,int nitem,int nout,real **c1,
 121                             real dt,unsigned long mode,int nrestart,
 122                             bool bAver,bool bFour,bool bNormalize,
 123                             bool bVerbose,real tbeginfit,real tendfit,
 124                             int nfitparm,int nskip);
 125 /*
 126  * do_autocorr calculates autocorrelation functions for many things.
 127  * It takes a 2 d array containing nitem arrays of length nframes
 128  * for each item the ACF is calculated.
 129  *
 130  * A number of "modes" exist for computation of the ACF
 131  *
 132  * if (mode == eacNormal) {
 133  *   C(t) = < X (tau) * X (tau+t) >
 134  * }
 135  * else if (mode == eacCos) {
 136  *   C(t) = < cos (X(tau) - X(tau+t)) >
 137  * }
 138  * else if (mode == eacVector) {
 139  *   C(t) = < X(tau) * X(tau+t)
 140  * }
 141  * else if (mode == eacP1) {
 142  *   C(t) = < cos (X(tau) * X(tau+t) >
 143  * }
 144  * else if (mode == eacP2) {
 145  *   C(t) = 1/2 * < 3 cos (X(tau) * X(tau+t) - 1 >
 146  * }
 147  * else if (mode == eacRcross) {
 148  *   C(t) = < ( X(tau) * X(tau+t) )^2 >
 149  * }
 150  *
 151  * For modes eacVector, eacP1, eacP2 and eacRcross the input should be
 152  * 3 x nframes long, where each triplet is taken as a 3D vector
 153  *
 154  * For mode eacCos inputdata must be in radians, not degrees!
 155  *
 156  * Other parameters are:
 157  *
 158  * fn is output filename (.xvg) where the correlation function(s) are printed
 159  * title is the title in the output file
 160  * nframes is the number of frames in the time series
 161  * nitem is the number of items
 162  * c1           is an array of dimension [ 0 .. nitem-1 ] [ 0 .. nframes-1 ]
 163  *              on output, this array is filled with the correlation function
 164  *              to reduce storage
 165  * nrestart     is the number of steps between restarts for direct ACFs
 166  *              (i.e. without FFT) When set to 1 all points are used as
 167  *              time origin for averaging
 168  * dt           is the time between frames
 169  * bAver        If set, all ndih C(t) functions are averaged into a single
 170  *              C(t)
 171  * bFour        If set, will use fast fourier transform (FFT) for evaluating
 172  *              the ACF
 173  * bNormalize   If set, all ACFs will be normalized to start at 0
 174  * nskip        Determines whether steps a re skipped in the output
 175  */
 176
 177 typedef struct {
 178   char *name;   /* Description of the J coupling constant */
 179   real A,B,C;   /* Karplus coefficients */
 180   real offset;  /* Offset for dihedral angle in histogram (e.g. -M_PI/3) */
 181   real Jc;      /* Resulting Jcoupling */
 182   real Jcsig;   /* Standard deviation in Jc */
 183 } t_karplus;
 184
 185 extern void calc_distribution_props(int nh,int histo[],
 186                                     real start,int  nkkk, t_karplus kkk[],
 187                                     real *S2);
 188 /* This routine takes a dihedral distribution and calculates
 189  * coupling constants and dihedral order parameters of it.
 190  *
 191  * nh      is the number of points
 192  * histo   is the array of datapoints which is assumed to span
 193  *         2 M_PI radians
 194  * start   is the starting angle of the histogram, this can be either 0
 195  *         or -M_PI
 196  * nkkk    is the number of karplus sets (multiple coupling constants may be
 197  *         derived from a single angle)
 198  * kkk     are the constants for calculating J coupling constants using a
 199  *         Karplus equation according to
 200  *
 201  *                  2
 202  *         J = A cos theta + B cos theta + C
 203  *
 204  *         where theta is phi - offset (phi is the angle in the histogram)
 205  * offset  is subtracted from phi before substitution in the Karplus
 206  *         equation
 207  * S2      is the resulting dihedral order parameter
 208  *
 209  */
 210
 211
 212 /***********************************************
 213  *
 214  *     F I T   R O U T I N E S
 215  *
 216  ***********************************************/
 217 extern void do_expfit(int ndata,real c1[],real dt,
 218                       real begintimefit,real endtimefit);
 219
 220 extern void expfit(int n, real x[], real y[], real Dy[],
 221                    real *a, real *sa,
 222                    real *b, real *sb);
 223 /* This procedure fits y=exp(a+bx) for n (x,y) pairs to determine a and b.
 224  * The uncertainties in the y values must be in the vector Dy.
 225  * The standard deviations of a and b, sa and sb, are also calculated.
 226  *
 227  * Routine from Computers in physics, 7(3) (1993), p. 280-285.
 228  */
 229
 230 extern void ana_dih_trans(char *fn_trans,char *fn_histo,
 231                           real **dih,int nframes,int nangles,
 232                           char *grpname,real t0,real dt,bool bRb);
 233 /*
 234  * Analyse dihedral transitions, by counting transitions per dihedral
 235  * and per frame. The total number of transitions is printed to
 236  * stderr, as well as the average time between transitions.
 237  *
 238  * is wrapper to low_ana_dih_trans, which also passes in and out the
 239      number of transitions per dihedral per residue. that uses struc dlist
 240      which is not external, so pp2shift.h must be included.
 241
 242  * Dihedrals are supposed to be in either of three minima,
 243  * (trans, gauche+, gauche-)
 244  *
 245  * fn_trans  output file name for #transitions per timeframe
 246  * fn_histo  output file name for transition time histogram
 247  * dih       the actual dihedral angles
 248  * nframes   number of times frames
 249  * nangles   number of angles
 250  * grpname   a string for the header of plots
 251  * t0,dt     starting time resp. time between time frames
 252  * bRb       determines whether the polymer convention is used
 253  *           (trans = 0)
 254  */
 255
 256 extern void low_ana_dih_trans(bool bTrans, char *fn_trans,
 257                               bool bHisto, char *fn_histo, int maxchi,
 258                               real **dih, int nlist, t_dlist dlist[], int nframes,
 259                               int nangles, char *grpname, int xity[],
 260                               real t0, real dt, bool bRb, real core_frac);
 261   /* as above but passes dlist so can copy occupancies into it, and xity[]
 262    *  (1..nangles, corresp to dih[this][], so can have non-3 multiplicity of
 263    * rotamers. Also production of xvg output files is conditional
 264    * and the fractional width of each rotamer can be set ie for a 3 fold
 265    * dihedral with core_frac = 0.5 only the central 60 degrees is assigned
 266    * to each rotamer, the rest goes to rotamer zero */
 267
 268
 269
 270 extern void read_ang_dih(char *trj_fn,char *tpb_fn,
 271                          bool bAngles,bool bSaveAll,bool bRb,
 272                          int maxangstat,int angstat[],
 273                          int *nframes,real **time,
 274                          int isize,atom_id index[],
 275                          real **trans_frac,
 276                          real **aver_angle,
 277                          real *dih[]);
 278 /*
 279  * Read a trajectory and calculate angles and dihedrals.
 280  *
 281  * trj_fn      file name of trajectory
 282  * tpb_fn      file name of tpb file
 283  * bAngles     do we have to read angles or dihedrals
 284  * bSaveAll    do we have to store all in the dih array
 285  * bRb         do we have Ryckaert-Bellemans dihedrals (trans = 0)
 286  * maxangstat  number of entries in distribution array
 287  * angstat     angle distribution
 288  * *nframes    number of frames read
 289  * time        simulation time at each time frame
 290  * isize       number of entries in the index, when angles 3*number of angles
 291  *             else 4*number of angles
 292  * index       atom numbers that define the angles or dihedrals
 293  *             (i,j,k) resp (i,j,k,l)
 294  * trans_frac  number of dihedrals in trans
 295  * aver_angle  average angle at each time frame
 296  * dih         all angles at each time frame
 297  */
 298
 299 extern void make_histo(FILE *log,
 300                        int ndata,real data[],int npoints,int histo[],
 301                        real minx,real maxx);
 302 /*
 303  * Make a histogram from data. The min and max of the data array can
 304  * be determined (if minx == 0 and maxx == 0)
 305  * and the index in the histogram is computed from
 306  * ind = npoints/(max(data) - min(data))
 307  *
 308  * log       write error output to this file
 309  * ndata     number of points in data
 310  * data      data points
 311  * npoints   number of points in histogram
 312  * histo     histogram array. This is NOT set to zero, to allow you
 313  *           to add multiple histograms
 314  * minx      start of the histogram
 315  * maxx      end of the histogram
 316  *           if both are 0, these values are computed by the routine itself
 317  */
 318
 319 extern void normalize_histo(int npoints,int histo[],real dx,real normhisto[]);
 320 /*
 321  * Normalize a histogram so that the integral over the histo is 1
 322  *
 323  * npoints    number of points in the histo array
 324  * histo      input histogram
 325  * dx         distance between points on the X-axis
 326  * normhisto  normalized output histogram
 327  */
 328
 329 extern real fit_function(int eFitFn,real *parm,real x);
 330 /* Returns the value of fit function eFitFn at x */
 331
 332 /* Use Levenberg-Marquardt method to fit to a nfitparm parameter exponential */
 333 /* or to a transverse current autocorrelation function */
 334 /* Or: "There is no KILL like OVERKILL", Dr. Ir. D. van der Spoel */
 335 extern real do_lmfit(int ndata,real c1[],real sig[],real dt,real *x,
 336                      real begintimefit,real endtimefit,bool bVerbose,
 337                      int eFitFn,real fitparms[],int fix);
 338 /* Returns integral.
 339  * If x == NULL, the timestep dt will be used to create a time axis.
 340  * fix fixes fit parameter i at it's starting value, when the i'th bit
 341  * of fix is set.
 342  */
 343
 344 extern real print_and_integrate(FILE *fp,int n,real dt,
 345                                 real c[],real *fit,int nskip);
 346 /* Integrate the data in c[] from 0 to n using trapezium rule.
 347  * If fp != NULL output is written to it
 348  * nskip determines whether all elements are written to the output file
 349  * (written when i % nskip == 0)
 350  * If fit != NULL the fit is also written.
 351  */
 352
 353 extern int get_acfnout(void);
 354 /* Return the output length for the correlation function
 355  * Works only AFTER do_auto_corr has been called!
 356  */
 357
 358 extern int get_acffitfn(void);
 359 /* Return the fit function type.
 360  * Works only AFTER do_auto_corr has been called!
 361  */
 362
 363
 364 /* Least squares method of computing statistics */
 365 typedef struct {
 366   double yy,yx,xx,sx,sy;
 367   int    np;
 368 } t_lsq;
 369
 370 extern void init_lsq(t_lsq *lsq);
 371 extern void done_lsq(t_lsq *lsq);
 372 extern void add_lsq(t_lsq *lsq,real x,real y);
 373 extern void get_lsq_ab(t_lsq *lsq,real *a,real *b);
 374 extern int  npoints_lsq(t_lsq *lsq);
 375 extern real aver_lsq(t_lsq *lsq);
 376 extern real sigma_lsq(t_lsq *lsq);
 377 extern real error_lsq(t_lsq *lsq);
 378
 379   /* Routines from pp2shift (anadih.c etc.) */
 380
 381 extern void do_pp2shifts(FILE *fp,int nframes,
 382                          int nlist,t_dlist dlist[],real **dih);
 383
 384 extern bool has_dihedral(int Dih,t_dlist *dl);
 385
 386 extern t_dlist *mk_dlist(FILE *log,
 387                          t_atoms *atoms, int *nlist,
 388                          bool bPhi, bool bPsi, bool bChi, bool bHChi,
 389                          int maxchi,int r0,int naa,char **aa);
 390
 391 extern void pr_dlist(FILE *fp,int nl,t_dlist dl[],real dt,  int printtype,
 392 bool bPhi, bool bPsi,bool bChi,bool bOmega, int maxchi);
 393
 394 extern int pr_trans(FILE *fp,int nl,t_dlist dl[],real dt,int Xi);
 395
 396 extern void mk_chi_lookup (int **lookup, int maxchi, real **dih,
 397                            int nlist, t_dlist dlist[]) ;
 398
 399 extern void mk_multiplicity_lookup (int *xity, int maxchi, real **dih,
 400                                     int nlist, t_dlist dlist[],int nangle) ;
 401
 402 extern void get_chi_product_traj (real **dih,int nframes,int nangles,
 403                                   int nlist,int maxchi, t_dlist dlist[], real time[],
 404                                   int **lookup,int *xity,bool bRb,bool bNormalize,
 405                                   real core_frac,bool bAll,char *fnall);
 406
 407 extern void print_one (char *base,char *name,char *title, char *ylabel,
 408                       int nf,real time[],real data[]);
 409
 410 #ifdef CPLUSPLUS
 411              }
 412 #endif
 413
 414 #endif
 415
 416