include/gstat.h

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  *                        VERSION 3.2.0
  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.
  14
  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.
  19  *
  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.
  26  *
  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.
  29  *
  30  * For more info, check our website at http://www.gromacs.org
  31  *
  32  * And Hey:
  33  * Gromacs Runs On Most of All Computer Systems
  34  */
  35
  36 #ifndef _gstat_h
  37 #define _gstat_h
  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 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 #define eacIden   (1<<9)
  66
  67 enum {
  68   effnNONE, effnEXP1, effnEXP2, effnEXP3,   effnVAC,
  69   effnEXP5, effnEXP7, effnEXP9, effnERF, effnERREST, effnNR
  70 };
  71
  72 /* must correspond with 'leg' g_chi.c:727 */
  73 enum { edPhi=0, edPsi, edOmega, edChi1, edChi2, edChi3, edChi4, edChi5, edChi6, edMax };
  74
  75 enum { edPrintST=0,edPrintRO } ;
  76
  77 #define NHISTO 360
  78 #define NONCHI 3
  79 #define MAXCHI edMax-NONCHI
  80 #define NROT 4  /* number of rotamers: 1=g(-), 2=t, 3=g(+), 0=other */
  81
  82 typedef struct {
  83   int minO,minC,H,N,C,O,Cn[MAXCHI+3];
  84 } t_dihatms; /* Cn[0]=N, Cn[1]=Ca, Cn[2]=Cb etc. */
  85
  86 typedef struct {
  87   char name[12];
  88   int  resnr;
  89   int  index;       /* Index for amino acids (histograms) */
  90   int  j0[edMax];   /* Index in dih array (phi angle is first...) */
  91   t_dihatms  atm;
  92   int  b[edMax];
  93   int  ntr[edMax];
  94   real S2[edMax];
  95   real rot_occ[edMax][NROT];
  96
  97 } t_dlist;
  98
  99 extern const int  nfp_ffn[effnNR];
 100
 101 extern const char *s_ffn[effnNR+2];
 102
 103 extern const char *longs_ffn[effnNR];
 104
 105 int sffn2effn(const char **sffn);
 106 /* Returns the ffn enum corresponding to the selected enum option in sffn */
 107
 108 t_pargs *add_acf_pargs(int *npargs,t_pargs *pa);
 109 /* Add options for autocorr to the current set of options.
 110  * *npargs must be initialised to the number of elements in pa,
 111  * it will be incremented appropriately.
 112  */
 113
 114 void cross_corr(int n,real f[],real g[],real corr[]);
 115 /* Simple minded cross correlation algorithm */
 116
 117 real fit_acf(int ncorr,int fitfn,const output_env_t oenv,gmx_bool bVerbose,
 118                     real tbeginfit,real tendfit,real dt,real c1[],real *fit);
 119   /* Fit an ACF to a given function */
 120
 121 void do_autocorr(const char *fn,const output_env_t oenv,
 122                         const char *title,
 123                         int nframes,int nitem,real **c1,
 124                         real dt,unsigned long mode,gmx_bool bAver);
 125 /* Calls low_do_autocorr (see below). After calling add_acf_pargs */
 126
 127 void low_do_autocorr(const char *fn,const output_env_t oenv,
 128                             const char *title, int  nframes,int nitem,
 129                             int nout,real **c1, real dt,unsigned long mode,
 130                             int nrestart, gmx_bool bAver,gmx_bool bNormalize,
 131                             gmx_bool bVerbose,real tbeginfit,real tendfit,
 132                             int nfitparm,int nskip);
 133 /*
 134  * do_autocorr calculates autocorrelation functions for many things.
 135  * It takes a 2 d array containing nitem arrays of length nframes
 136  * for each item the ACF is calculated.
 137  *
 138  * A number of "modes" exist for computation of the ACF
 139  *
 140  * if (mode == eacNormal) {
 141  *   C(t) = < X (tau) * X (tau+t) >
 142  * }
 143  * else if (mode == eacCos) {
 144  *   C(t) = < cos (X(tau) - X(tau+t)) >
 145  * }
 146  * else if (mode == eacIden) { **not fully supported yet**
 147  *   C(t) = < (X(tau) == X(tau+t)) >
 148  * }
 149  * else if (mode == eacVector) {
 150  *   C(t) = < X(tau) * X(tau+t)
 151  * }
 152  * else if (mode == eacP1) {
 153  *   C(t) = < cos (X(tau) * X(tau+t) >
 154  * }
 155  * else if (mode == eacP2) {
 156  *   C(t) = 1/2 * < 3 cos (X(tau) * X(tau+t) - 1 >
 157  * }
 158  * else if (mode == eacRcross) {
 159  *   C(t) = < ( X(tau) * X(tau+t) )^2 >
 160  * }
 161  *
 162  * For modes eacVector, eacP1, eacP2 and eacRcross the input should be
 163  * 3 x nframes long, where each triplet is taken as a 3D vector
 164  *
 165  * For mode eacCos inputdata must be in radians, not degrees!
 166  *
 167  * Other parameters are:
 168  *
 169  * fn is output filename (.xvg) where the correlation function(s) are printed
 170  * title is the title in the output file
 171  * nframes is the number of frames in the time series
 172  * nitem is the number of items
 173  * c1           is an array of dimension [ 0 .. nitem-1 ] [ 0 .. nframes-1 ]
 174  *              on output, this array is filled with the correlation function
 175  *              to reduce storage
 176  * nrestart     is the number of steps between restarts for direct ACFs
 177  *              (i.e. without FFT) When set to 1 all points are used as
 178  *              time origin for averaging
 179  * dt           is the time between frames
 180  * bAver        If set, all ndih C(t) functions are averaged into a single
 181  *              C(t)
 182  * (bFour       If set, will use fast fourier transform (FFT) for evaluating
 183  *              the ACF: removed option, now on the command line only)
 184  * bNormalize   If set, all ACFs will be normalized to start at 0
 185  * nskip        Determines whether steps a re skipped in the output
 186  */
 187
 188 typedef struct {
 189   const char *name; /* Description of the J coupling constant */
 190   real A,B,C;   /* Karplus coefficients */
 191   real offset;  /* Offset for dihedral angle in histogram (e.g. -M_PI/3) */
 192   real Jc;      /* Resulting Jcoupling */
 193   real Jcsig;   /* Standard deviation in Jc */
 194 } t_karplus;
 195
 196 void calc_distribution_props(int nh,int histo[],
 197                                     real start,int  nkkk, t_karplus kkk[],
 198                                     real *S2);
 199 /* This routine takes a dihedral distribution and calculates
 200  * coupling constants and dihedral order parameters of it.
 201  *
 202  * nh      is the number of points
 203  * histo   is the array of datapoints which is assumed to span
 204  *         2 M_PI radians
 205  * start   is the starting angle of the histogram, this can be either 0
 206  *         or -M_PI
 207  * nkkk    is the number of karplus sets (multiple coupling constants may be
 208  *         derived from a single angle)
 209  * kkk     are the constants for calculating J coupling constants using a
 210  *         Karplus equation according to
 211  *
 212  *                  2
 213  *         J = A cos theta + B cos theta + C
 214  *
 215  *         where theta is phi - offset (phi is the angle in the histogram)
 216  * offset  is subtracted from phi before substitution in the Karplus
 217  *         equation
 218  * S2      is the resulting dihedral order parameter
 219  *
 220  */
 221
 222
 223 /***********************************************
 224  *
 225  *     F I T   R O U T I N E S
 226  *
 227  ***********************************************/
 228 void do_expfit(int ndata,real c1[],real dt,
 229                       real begintimefit,real endtimefit);
 230
 231 void expfit(int n, real x[], real y[], real Dy[],
 232                    real *a, real *sa,
 233                    real *b, real *sb);
 234 /* This procedure fits y=exp(a+bx) for n (x,y) pairs to determine a and b.
 235  * The uncertainties in the y values must be in the vector Dy.
 236  * The standard deviations of a and b, sa and sb, are also calculated.
 237  *
 238  * Routine from Computers in physics, 7(3) (1993), p. 280-285.
 239  */
 240
 241 void ana_dih_trans(const char *fn_trans,const char *fn_histo,
 242                           real **dih,int nframes,int nangles,
 243                           const char *grpname,real t0,real dt,gmx_bool bRb,
 244                           const output_env_t oenv);
 245 /*
 246  * Analyse dihedral transitions, by counting transitions per dihedral
 247  * and per frame. The total number of transitions is printed to
 248  * stderr, as well as the average time between transitions.
 249  *
 250  * is wrapper to low_ana_dih_trans, which also passes in and out the
 251      number of transitions per dihedral per residue. that uses struc dlist
 252      which is not external, so pp2shift.h must be included.
 253
 254  * Dihedrals are supposed to be in either of three minima,
 255  * (trans, gauche+, gauche-)
 256  *
 257  * fn_trans  output file name for #transitions per timeframe
 258  * fn_histo  output file name for transition time histogram
 259  * dih       the actual dihedral angles
 260  * nframes   number of times frames
 261  * nangles   number of angles
 262  * grpname   a string for the header of plots
 263  * t0,dt     starting time resp. time between time frames
 264  * bRb       determines whether the polymer convention is used
 265  *           (trans = 0)
 266  */
 267
 268 void low_ana_dih_trans(gmx_bool bTrans, const char *fn_trans,
 269                               gmx_bool bHisto, const char *fn_histo, int maxchi,
 270                               real **dih, int nlist, t_dlist dlist[],
 271                               int nframes, int nangles, const char *grpname,
 272                                int multiplicity[], real t0, real dt, gmx_bool bRb,
 273                               real core_frac, const output_env_t oenv);
 274   /* as above but passes dlist so can copy occupancies into it, and multiplicity[]
 275    *  (1..nangles, corresp to dih[this][], so can have non-3 multiplicity of
 276    * rotamers. Also production of xvg output files is conditional
 277    * and the fractional width of each rotamer can be set ie for a 3 fold
 278    * dihedral with core_frac = 0.5 only the central 60 degrees is assigned
 279    * to each rotamer, the rest goes to rotamer zero */
 280
 281
 282
 283 void read_ang_dih(const char *trj_fn,
 284                          gmx_bool bAngles,gmx_bool bSaveAll,gmx_bool bRb,gmx_bool bPBC,
 285                          int maxangstat,int angstat[],
 286                          int *nframes,real **time,
 287                          int isize,atom_id index[],
 288                          real **trans_frac,
 289                          real **aver_angle,
 290                          real *dih[],
 291                          const output_env_t oenv);
 292 /*
 293  * Read a trajectory and calculate angles and dihedrals.
 294  *
 295  * trj_fn      file name of trajectory
 296  * tpb_fn      file name of tpb file
 297  * bAngles     do we have to read angles or dihedrals
 298  * bSaveAll    do we have to store all in the dih array
 299  * bRb         do we have Ryckaert-Bellemans dihedrals (trans = 0)
 300  * bPBC        compute angles module 2 Pi
 301  * maxangstat  number of entries in distribution array
 302  * angstat     angle distribution
 303  * *nframes    number of frames read
 304  * time        simulation time at each time frame
 305  * isize       number of entries in the index, when angles 3*number of angles
 306  *             else 4*number of angles
 307  * index       atom numbers that define the angles or dihedrals
 308  *             (i,j,k) resp (i,j,k,l)
 309  * trans_frac  number of dihedrals in trans
 310  * aver_angle  average angle at each time frame
 311  * dih         all angles at each time frame
 312  */
 313
 314 void make_histo(FILE *log,
 315                        int ndata,real data[],int npoints,int histo[],
 316                        real minx,real maxx);
 317 /*
 318  * Make a histogram from data. The min and max of the data array can
 319  * be determined (if minx == 0 and maxx == 0)
 320  * and the index in the histogram is computed from
 321  * ind = npoints/(max(data) - min(data))
 322  *
 323  * log       write error output to this file
 324  * ndata     number of points in data
 325  * data      data points
 326  * npoints   number of points in histogram
 327  * histo     histogram array. This is NOT set to zero, to allow you
 328  *           to add multiple histograms
 329  * minx      start of the histogram
 330  * maxx      end of the histogram
 331  *           if both are 0, these values are computed by the routine itself
 332  */
 333
 334 void normalize_histo(int npoints,int histo[],real dx,real normhisto[]);
 335 /*
 336  * Normalize a histogram so that the integral over the histo is 1
 337  *
 338  * npoints    number of points in the histo array
 339  * histo      input histogram
 340  * dx         distance between points on the X-axis
 341  * normhisto  normalized output histogram
 342  */
 343
 344 real fit_function(int eFitFn,real *parm,real x);
 345 /* Returns the value of fit function eFitFn at x */
 346
 347 /* Use Levenberg-Marquardt method to fit to a nfitparm parameter exponential */
 348 /* or to a transverse current autocorrelation function */
 349 /* Or: "There is no KILL like OVERKILL", Dr. Ir. D. van der Spoel */
 350 real do_lmfit(int ndata,real c1[],real sig[],real dt,real *x,
 351                      real begintimefit,real endtimefit,const output_env_t oenv,
 352                      gmx_bool bVerbose, int eFitFn,real fitparms[],int fix);
 353 /* Returns integral.
 354  * If x == NULL, the timestep dt will be used to create a time axis.
 355  * fix fixes fit parameter i at it's starting value, when the i'th bit
 356  * of fix is set.
 357  */
 358
 359 real evaluate_integral(int n,real x[],real y[],real dy[],
 360                               real aver_start,real *stddev);
 361 /* Integrate data in y, and, if given, use dy as weighting
 362  * aver_start should be set to a value where the function has
 363  * converged to 0.
 364  */
 365
 366 real print_and_integrate(FILE *fp,int n,real dt,
 367                                 real c[],real *fit,int nskip);
 368 /* Integrate the data in c[] from 0 to n using trapezium rule.
 369  * If fp != NULL output is written to it
 370  * nskip determines whether all elements are written to the output file
 371  * (written when i % nskip == 0)
 372  * If fit != NULL the fit is also written.
 373  */
 374
 375 int get_acfnout(void);
 376 /* Return the output length for the correlation function
 377  * Works only AFTER do_auto_corr has been called!
 378  */
 379
 380 int get_acffitfn(void);
 381 /* Return the fit function type.
 382  * Works only AFTER do_auto_corr has been called!
 383  */
 384
 385   /* Routines from pp2shift (anadih.c etc.) */
 386
 387 void do_pp2shifts(FILE *fp,int nframes,
 388                          int nlist,t_dlist dlist[],real **dih);
 389
 390 gmx_bool has_dihedral(int Dih,t_dlist *dl);
 391
 392 t_dlist *mk_dlist(FILE *log,
 393                          t_atoms *atoms, int *nlist,
 394                          gmx_bool bPhi, gmx_bool bPsi, gmx_bool bChi, gmx_bool bHChi,
 395                          int maxchi, int r0, gmx_residuetype_t rt);
 396
 397 void pr_dlist(FILE *fp,int nl,t_dlist dl[],real dt,  int printtype,
 398                      gmx_bool bPhi, gmx_bool bPsi,gmx_bool bChi,gmx_bool bOmega, int maxchi);
 399
 400 int pr_trans(FILE *fp,int nl,t_dlist dl[],real dt,int Xi);
 401
 402 void mk_chi_lookup (int **lookup, int maxchi, real **dih,
 403                            int nlist, t_dlist dlist[]) ;
 404
 405 void mk_multiplicity_lookup (int *multiplicity, int maxchi, real **dih,
 406                                     int nlist, t_dlist dlist[],int nangle) ;
 407
 408 void get_chi_product_traj (real **dih,int nframes,int nangles,
 409                                   int nlist,int maxchi, t_dlist dlist[],
 410                                   real time[], int **lookup,int *multiplicity,
 411                                   gmx_bool bRb,gmx_bool bNormalize,
 412                                   real core_frac,gmx_bool bAll,const char *fnall,
 413                                   const output_env_t oenv);
 414
 415 void print_one (const output_env_t oenv, const char *base,
 416                        const char *name,
 417                        const char *title, const char *ylabel,int nf,
 418                        real time[],real data[]);
 419
 420 /* Routines from g_hbond */
 421 void analyse_corr(int n,real t[],real ct[],real nt[],real kt[],
 422                          real sigma_ct[],real sigma_nt[],real sigma_kt[],
 423                          real fit_start,real temp,real smooth_tail_start,
 424                          const output_env_t oenv);
 425
 426 void compute_derivative(int nn,real x[],real y[],real dydx[]);
 427
 428 #ifdef __cplusplus
 429              }
 430 #endif
 431
 432 #endif
 433
 434