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,2021, 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 struct t_dlist
  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
  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                        gmx::ArrayRef<t_dlist>  dlist,
 163                        int                     nframes,
 164                        int                     nangles,
 165                        const char*             grpname,
 166                        int                     multiplicity[],
 167                        real*                   time,
 168                        gmx_bool                bRb,
 169                        real                    core_frac,
 170                        const gmx_output_env_t* oenv);
 171 /* as above but passes dlist so can copy occupancies into it, and multiplicity[]
 172  *  (1..nangles, corresp to dih[this][], so can have non-3 multiplicity of
 173  * rotamers. Also production of xvg output files is conditional
 174  * and the fractional width of each rotamer can be set ie for a 3 fold
 175  * dihedral with core_frac = 0.5 only the central 60 degrees is assigned
 176  * to each rotamer, the rest goes to rotamer zero */
 177
 178
 179 void read_ang_dih(const char*             trj_fn,
 180                   gmx_bool                bAngles,
 181                   gmx_bool                bSaveAll,
 182                   gmx_bool                bRb,
 183                   gmx_bool                bPBC,
 184                   int                     maxangstat,
 185                   int                     angstat[],
 186                   int*                    nframes,
 187                   real**                  time,
 188                   int                     isize,
 189                   int                     index[],
 190                   real**                  trans_frac,
 191                   real**                  aver_angle,
 192                   real*                   dih[],
 193                   const gmx_output_env_t* oenv);
 194 /*
 195  * Read a trajectory and calculate angles and dihedrals.
 196  *
 197  * trj_fn      file name of trajectory
 198  * bAngles     do we have to read angles or dihedrals
 199  * bSaveAll    do we have to store all in the dih array
 200  * bRb         do we have Ryckaert-Bellemans dihedrals (trans = 0)
 201  * bPBC        compute angles module 2 Pi
 202  * maxangstat  number of entries in distribution array
 203  * angstat     angle distribution
 204  * *nframes    number of frames read
 205  * time        simulation time at each time frame
 206  * isize       number of entries in the index, when angles 3*number of angles
 207  *             else 4*number of angles
 208  * index       atom numbers that define the angles or dihedrals
 209  *             (i,j,k) resp (i,j,k,l)
 210  * trans_frac  number of dihedrals in trans
 211  * aver_angle  average angle at each time frame
 212  * dih         all angles at each time frame
 213  */
 214
 215 void make_histo(FILE* log, int ndata, real data[], int npoints, int histo[], real minx, real maxx);
 216 /*
 217  * Make a histogram from data. The min and max of the data array can
 218  * be determined (if minx == 0 and maxx == 0)
 219  * and the index in the histogram is computed from
 220  * ind = npoints/(max(data) - min(data))
 221  *
 222  * log       write error output to this file
 223  * ndata     number of points in data
 224  * data      data points
 225  * npoints   number of points in histogram
 226  * histo     histogram array. This is NOT set to zero, to allow you
 227  *           to add multiple histograms
 228  * minx      start of the histogram
 229  * maxx      end of the histogram
 230  *           if both are 0, these values are computed by the routine itself
 231  */
 232
 233 void normalize_histo(gmx::ArrayRef<const int> histo, real dx, gmx::ArrayRef<real> normhisto);
 234 /*
 235  * Normalize a histogram so that the integral over the histo is 1
 236  *
 237  * histo      input histogram
 238  * dx         distance between points on the X-axis
 239  * normhisto  normalized output histogram
 240  */
 241
 242 /* Routines from pp2shift (anadih.c etc.) */
 243
 244 void do_pp2shifts(FILE* fp, int nframes, gmx::ArrayRef<const t_dlist> dlist, real** dih);
 245
 246 gmx_bool has_dihedral(int Dih, const t_dlist& dlist);
 247
 248 std::vector<t_dlist> mk_dlist(FILE*          log,
 249                               const t_atoms* atoms,
 250                               gmx_bool       bPhi,
 251                               gmx_bool       bPsi,
 252                               gmx_bool       bChi,
 253                               gmx_bool       bHChi,
 254                               int            maxchi,
 255                               int            r0,
 256                               ResidueType*   rt);
 257
 258 void pr_dlist(FILE*                        fp,
 259               gmx::ArrayRef<const t_dlist> dlist,
 260               real                         dt,
 261               int                          printtype,
 262               gmx_bool                     bPhi,
 263               gmx_bool                     bPsi,
 264               gmx_bool                     bChi,
 265               gmx_bool                     bOmega,
 266               int                          maxchi);
 267
 268 void mk_chi_lookup(int** lookup, int maxchi, gmx::ArrayRef<const t_dlist> dlist);
 269
 270 void mk_multiplicity_lookup(int* multiplicity, int maxchi, gmx::ArrayRef<const t_dlist> dlist, int nangle);
 271
 272 void get_chi_product_traj(real**                       dih,
 273                           int                          nframes,
 274                           int                          maxchi,
 275                           gmx::ArrayRef<const t_dlist> dlist,
 276                           real                         time[],
 277                           int**                        lookup,
 278                           int*                         multiplicity,
 279                           gmx_bool                     bRb,
 280                           gmx_bool                     bNormalize,
 281                           real                         core_frac,
 282                           gmx_bool                     bAll,
 283                           const char*                  fnall,
 284                           const gmx_output_env_t*      oenv);
 285
 286 void print_one(const gmx_output_env_t* oenv,
 287                const char*             base,
 288                const char*             name,
 289                const char*             title,
 290                const char*             ylabel,
 291                int                     nf,
 292                real                    time[],
 293                real                    data[]);
 294
 295 /* Routines from g_hbond */
 296 void analyse_corr(int  n,
 297                   real t[],
 298                   real ct[],
 299                   real nt[],
 300                   real kt[],
 301                   real sigma_ct[],
 302                   real sigma_nt[],
 303                   real sigma_kt[],
 304                   real fit_start,
 305                   real temp);
 306
 307 void compute_derivative(int nn, const real x[], const real y[], real dydx[]);
 308
 309 #endif