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
  48 /* must correspond with 'leg' g_chi.c:727 */
  49 enum
  50 {
  51     edPhi = 0,
  52     edPsi,
  53     edOmega,
  54     edChi1,
  55     edChi2,
  56     edChi3,
  57     edChi4,
  58     edChi5,
  59     edChi6,
  60     edMax
  61 };
  62
  63 enum
  64 {
  65     edPrintST = 0,
  66     edPrintRO
  67 };
  68
  69 #define NHISTO 360
  70 #define NONCHI 3
  71 #define MAXCHI (edMax - NONCHI)
  72 #define NROT 4 /* number of rotamers: 1=g(-), 2=t, 3=g(+), 0=other */
  73
  74 typedef struct
  75 {
  76     int minCalpha, minC, H, N, C, O, Cn[MAXCHI + 3];
  77 } t_dihatms; /* Cn[0]=N, Cn[1]=Ca, Cn[2]=Cb etc. */
  78
  79 struct t_dlist
  80 {
  81     char        name[12];
  82     int         resnr;
  83     std::string residueName;
  84     int         j0[edMax]; /* Index in dih array (phi angle is first...) */
  85     t_dihatms   atm;
  86     int         b[edMax];
  87     int         ntr[edMax];
  88     real        S2[edMax];
  89     real        rot_occ[edMax][NROT];
  90 };
  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                        gmx::ArrayRef<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(gmx::ArrayRef<const int> histo, real dx, gmx::ArrayRef<real> normhisto);
 235 /*
 236  * Normalize a histogram so that the integral over the histo is 1
 237  *
 238  * histo      input histogram
 239  * dx         distance between points on the X-axis
 240  * normhisto  normalized output histogram
 241  */
 242
 243 /* Routines from pp2shift (anadih.c etc.) */
 244
 245 void do_pp2shifts(FILE* fp, int nframes, gmx::ArrayRef<const t_dlist> dlist, real** dih);
 246
 247 gmx_bool has_dihedral(int Dih, const t_dlist& dlist);
 248
 249 /*! \brief Describe the dihedrals in the residues of the \c atoms
 250  * structure
 251  *
 252  * Return a vector with a t_dlist entry for each residue in \c
 253  * atoms. The entry for a residue contains its name, its index within
 254  * the residues, and a mapping from chemical peptide atom names to
 255  * atom indices based on the atom names. Many fields of t_dlist are
 256  * not yet filled. */
 257 std::vector<t_dlist> mk_dlist(FILE*          log,
 258                               const t_atoms* atoms,
 259                               gmx_bool       bPhi,
 260                               gmx_bool       bPsi,
 261                               gmx_bool       bChi,
 262                               gmx_bool       bHChi,
 263                               int            maxchi,
 264                               int            r0);
 265
 266 void pr_dlist(FILE*                        fp,
 267               gmx::ArrayRef<const t_dlist> dlist,
 268               real                         dt,
 269               int                          printtype,
 270               gmx_bool                     bPhi,
 271               gmx_bool                     bPsi,
 272               gmx_bool                     bChi,
 273               gmx_bool                     bOmega,
 274               int                          maxchi);
 275
 276 void mk_chi_lookup(int** lookup, int maxchi, gmx::ArrayRef<const t_dlist> dlist);
 277
 278 void mk_multiplicity_lookup(int* multiplicity, int maxchi, gmx::ArrayRef<const t_dlist> dlist, int nangle);
 279
 280 void get_chi_product_traj(real**                       dih,
 281                           int                          nframes,
 282                           int                          maxchi,
 283                           gmx::ArrayRef<const 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