include/gstat.h

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