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 <string>
  42
  43 #include "gromacs/commandline/pargs.h"
  44 #include "gromacs/topology/index.h"
  45
  46 struct gmx_output_env_t;
  47 class ResidueTypeMap;
  48
  49 /* must correspond with 'leg' g_chi.c:727 */
  50 enum
  51 {
  52     edPhi = 0,
  53     edPsi,
  54     edOmega,
  55     edChi1,
  56     edChi2,
  57     edChi3,
  58     edChi4,
  59     edChi5,
  60     edChi6,
  61     edMax
  62 };
  63
  64 enum
  65 {
  66     edPrintST = 0,
  67     edPrintRO
  68 };
  69
  70 #define NHISTO 360
  71 #define NONCHI 3
  72 #define MAXCHI (edMax - NONCHI)
  73 #define NROT 4 /* number of rotamers: 1=g(-), 2=t, 3=g(+), 0=other */
  74
  75 typedef struct
  76 {
  77     int minCalpha, minC, H, N, C, O, Cn[MAXCHI + 3];
  78 } t_dihatms; /* Cn[0]=N, Cn[1]=Ca, Cn[2]=Cb etc. */
  79
  80 struct t_dlist
  81 {
  82     char        name[12];
  83     int         resnr;
  84     std::string residueName;
  85     int         j0[edMax]; /* Index in dih array (phi angle is first...) */
  86     t_dihatms   atm;
  87     int         b[edMax];
  88     int         ntr[edMax];
  89     real        S2[edMax];
  90     real        rot_occ[edMax][NROT];
  91 };
  92
  93 typedef struct
  94 {
  95     const char* name;    /* Description of the J coupling constant */
  96     real        A, B, C; /* Karplus coefficients */
  97     real        offset;  /* Offset for dihedral angle in histogram (e.g. -M_PI/3) */
  98     real        Jc;      /* Resulting Jcoupling */
  99     real        Jcsig;   /* Standard deviation in Jc */
 100 } t_karplus;
 101
 102 void calc_distribution_props(int nh, const int histo[], real start, int nkkk, t_karplus kkk[], real* S2);
 103 /* This routine takes a dihedral distribution and calculates
 104  * coupling constants and dihedral order parameters of it.
 105  *
 106  * nh      is the number of points
 107  * histo   is the array of datapoints which is assumed to span
 108  *         2 M_PI radians
 109  * start   is the starting angle of the histogram, this can be either 0
 110  *         or -M_PI
 111  * nkkk    is the number of karplus sets (multiple coupling constants may be
 112  *         derived from a single angle)
 113  * kkk     are the constants for calculating J coupling constants using a
 114  *         Karplus equation according to
 115  *
 116  *                  2
 117  *         J = A cos theta + B cos theta + C
 118  *
 119  *         where theta is phi - offset (phi is the angle in the histogram)
 120  * offset  is subtracted from phi before substitution in the Karplus
 121  *         equation
 122  * S2      is the resulting dihedral order parameter
 123  *
 124  */
 125
 126 void ana_dih_trans(const char*             fn_trans,
 127                    const char*             fn_histo,
 128                    real**                  dih,
 129                    int                     nframes,
 130                    int                     nangles,
 131                    const char*             grpname,
 132                    real*                   time,
 133                    gmx_bool                bRb,
 134                    const gmx_output_env_t* oenv);
 135 /*
 136  * Analyse dihedral transitions, by counting transitions per dihedral
 137  * and per frame. The total number of transitions is printed to
 138  * stderr, as well as the average time between transitions.
 139  *
 140  * is wrapper to low_ana_dih_trans, which also passes in and out the
 141      number of transitions per dihedral per residue. that uses struc dlist
 142      which is not external, so pp2shift.h must be included.
 143
 144  * Dihedrals are supposed to be in either of three minima,
 145  * (trans, gauche+, gauche-)
 146  *
 147  * fn_trans  output file name for #transitions per timeframe
 148  * fn_histo  output file name for transition time histogram
 149  * dih       the actual dihedral angles
 150  * nframes   number of times frames
 151  * nangles   number of angles
 152  * grpname   a string for the header of plots
 153  * time      array (size nframes) of times of trajectory frames
 154  * bRb       determines whether the polymer convention is used
 155  *           (trans = 0)
 156  */
 157
 158 void low_ana_dih_trans(gmx_bool                bTrans,
 159                        const char*             fn_trans,
 160                        gmx_bool                bHisto,
 161                        const char*             fn_histo,
 162                        int                     maxchi,
 163                        real**                  dih,
 164                        gmx::ArrayRef<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(gmx::ArrayRef<const int> histo, real dx, gmx::ArrayRef<real> normhisto);
 236 /*
 237  * Normalize a histogram so that the integral over the histo is 1
 238  *
 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, gmx::ArrayRef<const t_dlist> dlist, real** dih);
 247
 248 gmx_bool has_dihedral(int Dih, const t_dlist& dlist);
 249
 250 std::vector<t_dlist> mk_dlist(FILE*           log,
 251                               const t_atoms*  atoms,
 252                               gmx_bool        bPhi,
 253                               gmx_bool        bPsi,
 254                               gmx_bool        bChi,
 255                               gmx_bool        bHChi,
 256                               int             maxchi,
 257                               int             r0,
 258                               ResidueTypeMap* rt);
 259
 260 void pr_dlist(FILE*                        fp,
 261               gmx::ArrayRef<const t_dlist> dlist,
 262               real                         dt,
 263               int                          printtype,
 264               gmx_bool                     bPhi,
 265               gmx_bool                     bPsi,
 266               gmx_bool                     bChi,
 267               gmx_bool                     bOmega,
 268               int                          maxchi);
 269
 270 void mk_chi_lookup(int** lookup, int maxchi, gmx::ArrayRef<const t_dlist> dlist);
 271
 272 void mk_multiplicity_lookup(int* multiplicity, int maxchi, gmx::ArrayRef<const t_dlist> dlist, int nangle);
 273
 274 void get_chi_product_traj(real**                       dih,
 275                           int                          nframes,
 276                           int                          maxchi,
 277                           gmx::ArrayRef<const t_dlist> dlist,
 278                           real                         time[],
 279                           int**                        lookup,
 280                           int*                         multiplicity,
 281                           gmx_bool                     bRb,
 282                           gmx_bool                     bNormalize,
 283                           real                         core_frac,
 284                           gmx_bool                     bAll,
 285                           const char*                  fnall,
 286                           const gmx_output_env_t*      oenv);
 287
 288 void print_one(const gmx_output_env_t* oenv,
 289                const char*             base,
 290                const char*             name,
 291                const char*             title,
 292                const char*             ylabel,
 293                int                     nf,
 294                real                    time[],
 295                real                    data[]);
 296
 297 /* Routines from g_hbond */
 298 void analyse_corr(int  n,
 299                   real t[],
 300                   real ct[],
 301                   real nt[],
 302                   real kt[],
 303                   real sigma_ct[],
 304                   real sigma_nt[],
 305                   real sigma_kt[],
 306                   real fit_start,
 307                   real temp);
 308
 309 void compute_derivative(int nn, const real x[], const real y[], real dydx[]);
 310
 311 #endif