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