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