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