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,2015,2016,2018, 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 "gromacs/commandline/pargs.h"
  41 #include "gromacs/topology/index.h"
  42
  43 #ifdef __cplusplus
  44 extern "C" {
  45 #endif
  46
  47 struct gmx_output_env_t;
  48 struct gmx_residuetype_t;
  49
  50 /* must correspond with 'leg' g_chi.c:727 */
  51 enum {
  52     edPhi = 0, edPsi, edOmega, edChi1, edChi2, edChi3, edChi4, edChi5, edChi6, edMax
  53 };
  54
  55 enum {
  56     edPrintST = 0, edPrintRO
  57 };
  58
  59 #define NHISTO 360
  60 #define NONCHI 3
  61 #define MAXCHI (edMax-NONCHI)
  62 #define NROT 4  /* number of rotamers: 1=g(-), 2=t, 3=g(+), 0=other */
  63
  64 typedef struct {
  65     int minCalpha, minC, H, N, C, O, Cn[MAXCHI+3];
  66 } t_dihatms; /* Cn[0]=N, Cn[1]=Ca, Cn[2]=Cb etc. */
  67
  68 typedef struct {
  69     char       name[12];
  70     int        resnr;
  71     int        index;     /* Index for amino acids (histograms) */
  72     int        j0[edMax]; /* Index in dih array (phi angle is first...) */
  73     t_dihatms  atm;
  74     int        b[edMax];
  75     int        ntr[edMax];
  76     real       S2[edMax];
  77     real       rot_occ[edMax][NROT];
  78
  79 } t_dlist;
  80
  81 typedef struct {
  82     const char *name;    /* Description of the J coupling constant */
  83     real        A, B, C; /* Karplus coefficients */
  84     real        offset;  /* Offset for dihedral angle in histogram (e.g. -M_PI/3) */
  85     real        Jc;      /* Resulting Jcoupling */
  86     real        Jcsig;   /* Standard deviation in Jc */
  87 } t_karplus;
  88
  89 void calc_distribution_props(int nh, const int histo[],
  90                              real start, int  nkkk, t_karplus kkk[],
  91                              real *S2);
  92 /* This routine takes a dihedral distribution and calculates
  93  * coupling constants and dihedral order parameters of it.
  94  *
  95  * nh      is the number of points
  96  * histo   is the array of datapoints which is assumed to span
  97  *         2 M_PI radians
  98  * start   is the starting angle of the histogram, this can be either 0
  99  *         or -M_PI
 100  * nkkk    is the number of karplus sets (multiple coupling constants may be
 101  *         derived from a single angle)
 102  * kkk     are the constants for calculating J coupling constants using a
 103  *         Karplus equation according to
 104  *
 105  *                  2
 106  *         J = A cos theta + B cos theta + C
 107  *
 108  *         where theta is phi - offset (phi is the angle in the histogram)
 109  * offset  is subtracted from phi before substitution in the Karplus
 110  *         equation
 111  * S2      is the resulting dihedral order parameter
 112  *
 113  */
 114
 115 void ana_dih_trans(const char *fn_trans, const char *fn_histo,
 116                    real **dih, int nframes, int nangles,
 117                    const char *grpname, real *time, gmx_bool bRb,
 118                    const gmx_output_env_t *oenv);
 119 /*
 120  * Analyse dihedral transitions, by counting transitions per dihedral
 121  * and per frame. The total number of transitions is printed to
 122  * stderr, as well as the average time between transitions.
 123  *
 124  * is wrapper to low_ana_dih_trans, which also passes in and out the
 125      number of transitions per dihedral per residue. that uses struc dlist
 126      which is not external, so pp2shift.h must be included.
 127
 128  * Dihedrals are supposed to be in either of three minima,
 129  * (trans, gauche+, gauche-)
 130  *
 131  * fn_trans  output file name for #transitions per timeframe
 132  * fn_histo  output file name for transition time histogram
 133  * dih       the actual dihedral angles
 134  * nframes   number of times frames
 135  * nangles   number of angles
 136  * grpname   a string for the header of plots
 137  * time      array (size nframes) of times of trajectory frames
 138  * bRb       determines whether the polymer convention is used
 139  *           (trans = 0)
 140  */
 141
 142 void low_ana_dih_trans(gmx_bool bTrans, const char *fn_trans,
 143                        gmx_bool bHisto, const char *fn_histo, int maxchi,
 144                        real **dih, int nlist, t_dlist dlist[],
 145                        int nframes, int nangles, const char *grpname,
 146                        int multiplicity[], real *time, gmx_bool bRb,
 147                        real core_frac, const gmx_output_env_t *oenv);
 148 /* as above but passes dlist so can copy occupancies into it, and multiplicity[]
 149  *  (1..nangles, corresp to dih[this][], so can have non-3 multiplicity of
 150  * rotamers. Also production of xvg output files is conditional
 151  * and the fractional width of each rotamer can be set ie for a 3 fold
 152  * dihedral with core_frac = 0.5 only the central 60 degrees is assigned
 153  * to each rotamer, the rest goes to rotamer zero */
 154
 155
 156
 157 void read_ang_dih(const char *trj_fn,
 158                   gmx_bool bAngles, gmx_bool bSaveAll, gmx_bool bRb, gmx_bool bPBC,
 159                   int maxangstat, int angstat[],
 160                   int *nframes, real **time,
 161                   int isize, int index[],
 162                   real **trans_frac,
 163                   real **aver_angle,
 164                   real *dih[],
 165                   const gmx_output_env_t *oenv);
 166 /*
 167  * Read a trajectory and calculate angles and dihedrals.
 168  *
 169  * trj_fn      file name of trajectory
 170  * bAngles     do we have to read angles or dihedrals
 171  * bSaveAll    do we have to store all in the dih array
 172  * bRb         do we have Ryckaert-Bellemans dihedrals (trans = 0)
 173  * bPBC        compute angles module 2 Pi
 174  * maxangstat  number of entries in distribution array
 175  * angstat     angle distribution
 176  * *nframes    number of frames read
 177  * time        simulation time at each time frame
 178  * isize       number of entries in the index, when angles 3*number of angles
 179  *             else 4*number of angles
 180  * index       atom numbers that define the angles or dihedrals
 181  *             (i,j,k) resp (i,j,k,l)
 182  * trans_frac  number of dihedrals in trans
 183  * aver_angle  average angle at each time frame
 184  * dih         all angles at each time frame
 185  */
 186
 187 void make_histo(FILE *log,
 188                 int ndata, real data[], int npoints, int histo[],
 189                 real minx, real maxx);
 190 /*
 191  * Make a histogram from data. The min and max of the data array can
 192  * be determined (if minx == 0 and maxx == 0)
 193  * and the index in the histogram is computed from
 194  * ind = npoints/(max(data) - min(data))
 195  *
 196  * log       write error output to this file
 197  * ndata     number of points in data
 198  * data      data points
 199  * npoints   number of points in histogram
 200  * histo     histogram array. This is NOT set to zero, to allow you
 201  *           to add multiple histograms
 202  * minx      start of the histogram
 203  * maxx      end of the histogram
 204  *           if both are 0, these values are computed by the routine itself
 205  */
 206
 207 void normalize_histo(int npoints, const int histo[], real dx, real normhisto[]);
 208 /*
 209  * Normalize a histogram so that the integral over the histo is 1
 210  *
 211  * npoints    number of points in the histo array
 212  * histo      input histogram
 213  * dx         distance between points on the X-axis
 214  * normhisto  normalized output histogram
 215  */
 216
 217 /* Routines from pp2shift (anadih.c etc.) */
 218
 219 void do_pp2shifts(FILE *fp, int nframes,
 220                   int nlist, t_dlist dlist[], real **dih);
 221
 222 gmx_bool has_dihedral(int Dih, t_dlist *dl);
 223
 224 t_dlist *mk_dlist(FILE *log,
 225                   const t_atoms *atoms, int *nlist,
 226                   gmx_bool bPhi, gmx_bool bPsi, gmx_bool bChi, gmx_bool bHChi,
 227                   int maxchi, int r0, struct gmx_residuetype_t *rt);
 228
 229 void pr_dlist(FILE *fp, int nl, t_dlist dl[], real dt,  int printtype,
 230               gmx_bool bPhi, gmx_bool bPsi, gmx_bool bChi, gmx_bool bOmega, int maxchi);
 231
 232 int pr_trans(FILE *fp, int nl, t_dlist dl[], real dt, int Xi);
 233
 234 void mk_chi_lookup (int **lookup, int maxchi,
 235                     int nlist, t_dlist dlist[]);
 236
 237 void mk_multiplicity_lookup (int *multiplicity, int maxchi,
 238                              int nlist, t_dlist dlist[], int nangle);
 239
 240 void get_chi_product_traj (real **dih, int nframes,
 241                            int nlist, int maxchi, t_dlist dlist[],
 242                            real time[], int **lookup, int *multiplicity,
 243                            gmx_bool bRb, gmx_bool bNormalize,
 244                            real core_frac, gmx_bool bAll, const char *fnall,
 245                            const gmx_output_env_t *oenv);
 246
 247 void print_one (const gmx_output_env_t *oenv, const char *base,
 248                 const char *name,
 249                 const char *title, const char *ylabel, int nf,
 250                 real time[], real data[]);
 251
 252 /* Routines from g_hbond */
 253 void analyse_corr(int n, real t[], real ct[], real nt[], real kt[],
 254                   real sigma_ct[], real sigma_nt[], real sigma_kt[],
 255                   real fit_start, real temp);
 256
 257 void compute_derivative(int nn, const real x[], const real y[], real dydx[]);
 258
 259 #ifdef __cplusplus
 260 }
 261 #endif
 262
 263 #endif