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