src/gromacs/gmxana/gmx_angle.c

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  *                        VERSION 3.2.0
  10  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  11  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  12  * Copyright (c) 2001-2004, The GROMACS development team,
  13  * check out http://www.gromacs.org for more information.
  14
  15  * This program is free software; you can redistribute it and/or
  16  * modify it under the terms of the GNU General Public License
  17  * as published by the Free Software Foundation; either version 2
  18  * of the License, or (at your option) any later version.
  19  *
  20  * If you want to redistribute modifications, please consider that
  21  * scientific software is very special. Version control is crucial -
  22  * bugs must be traceable. We will be happy to consider code for
  23  * inclusion in the official distribution, but derived work must not
  24  * be called official GROMACS. Details are found in the README & COPYING
  25  * files - if they are missing, get the official version at www.gromacs.org.
  26  *
  27  * To help us fund GROMACS development, we humbly ask that you cite
  28  * the papers on the package - you can find them in the top README file.
  29  *
  30  * For more info, check our website at http://www.gromacs.org
  31  *
  32  * And Hey:
  33  * Green Red Orange Magenta Azure Cyan Skyblue
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38 #include <math.h>
  39 #include <string.h>
  40
  41 #include "sysstuff.h"
  42 #include "physics.h"
  43 #include "typedefs.h"
  44 #include "smalloc.h"
  45 #include "futil.h"
  46 #include "statutil.h"
  47 #include "copyrite.h"
  48 #include "vec.h"
  49 #include "index.h"
  50 #include "macros.h"
  51 #include "gmx_fatal.h"
  52 #include "xvgr.h"
  53 #include "gstat.h"
  54 #include "trnio.h"
  55 #include "gmx_ana.h"
  56
  57
  58 static void dump_dih_trn(int nframes, int nangles, real **dih, const char *fn,
  59                          real *time)
  60 {
  61     int       i, j, k, l, m, na;
  62     t_fileio *trn;
  63     rvec     *x;
  64     matrix    box = {{2, 0, 0}, {0, 2, 0}, {0, 0, 2}};
  65
  66     na = (nangles*2);
  67     if ((na % 3) != 0)
  68     {
  69         na = 1+na/3;
  70     }
  71     else
  72     {
  73         na = na/3;
  74     }
  75     printf("There are %d dihedrals. Will fill %d atom positions with cos/sin\n",
  76            nangles, na);
  77     snew(x, na);
  78     trn = open_trn(fn, "w");
  79     for (i = 0; (i < nframes); i++)
  80     {
  81         k = l = 0;
  82         for (j = 0; (j < nangles); j++)
  83         {
  84             for (m = 0; (m < 2); m++)
  85             {
  86                 x[k][l] = (m == 0) ? cos(dih[j][i]) : sin(dih[j][i]);
  87                 l++;
  88                 if (l == DIM)
  89                 {
  90                     l = 0;
  91                     k++;
  92                 }
  93             }
  94         }
  95         fwrite_trn(trn, i, time[i], 0, box, na, x, NULL, NULL);
  96     }
  97     close_trn(trn);
  98     sfree(x);
  99 }
 100
 101 int gmx_g_angle(int argc, char *argv[])
 102 {
 103     static const char *desc[] = {
 104         "[TT]g_angle[tt] computes the angle distribution for a number of angles",
 105         "or dihedrals.[PAR]",
 106         "With option [TT]-ov[tt], you can plot the average angle of",
 107         "a group of angles as a function of time. With the [TT]-all[tt] option,",
 108         "the first graph is the average and the rest are the individual angles.[PAR]",
 109         "With the [TT]-of[tt] option, [TT]g_angle[tt] also calculates the fraction of trans",
 110         "dihedrals (only for dihedrals) as function of time, but this is",
 111         "probably only fun for a select few.[PAR]",
 112         "With option [TT]-oc[tt], a dihedral correlation function is calculated.[PAR]",
 113         "It should be noted that the index file must contain",
 114         "atom triplets for angles or atom quadruplets for dihedrals.",
 115         "If this is not the case, the program will crash.[PAR]",
 116         "With option [TT]-or[tt], a trajectory file is dumped containing cos and",
 117         "sin of selected dihedral angles, which subsequently can be used as",
 118         "input for a principal components analysis using [TT]g_covar[tt].[PAR]",
 119         "Option [TT]-ot[tt] plots when transitions occur between",
 120         "dihedral rotamers of multiplicity 3 and [TT]-oh[tt]",
 121         "records a histogram of the times between such transitions,",
 122         "assuming the input trajectory frames are equally spaced in time."
 123     };
 124     static const char *opt[]    = { NULL, "angle", "dihedral", "improper", "ryckaert-bellemans", NULL };
 125     static gmx_bool    bALL     = FALSE, bChandler = FALSE, bAverCorr = FALSE, bPBC = TRUE;
 126     static real        binwidth = 1;
 127     t_pargs            pa[]     = {
 128         { "-type", FALSE, etENUM, {opt},
 129           "Type of angle to analyse" },
 130         { "-all",    FALSE,  etBOOL, {&bALL},
 131           "Plot all angles separately in the averages file, in the order of appearance in the index file." },
 132         { "-binwidth", FALSE, etREAL, {&binwidth},
 133           "binwidth (degrees) for calculating the distribution" },
 134         { "-periodic", FALSE, etBOOL, {&bPBC},
 135           "Print dihedral angles modulo 360 degrees" },
 136         { "-chandler", FALSE,  etBOOL, {&bChandler},
 137           "Use Chandler correlation function (N[trans] = 1, N[gauche] = 0) rather than cosine correlation function. Trans is defined as phi < -60 or phi > 60." },
 138         { "-avercorr", FALSE,  etBOOL, {&bAverCorr},
 139           "Average the correlation functions for the individual angles/dihedrals" }
 140     };
 141     static const char *bugs[] = {
 142         "Counting transitions only works for dihedrals with multiplicity 3"
 143     };
 144
 145     FILE              *out;
 146     real               tmp, dt;
 147     int                status, isize;
 148     atom_id           *index;
 149     char              *grpname;
 150     real               maxang, Jc, S2, norm_fac, maxstat;
 151     unsigned long      mode;
 152     int                nframes, maxangstat, mult, *angstat;
 153     int                i, j, total, nangles, natoms, nat2, first, last, angind;
 154     gmx_bool           bAver, bRb, bPeriodic,
 155                        bFrac,                           /* calculate fraction too?  */
 156                        bTrans,                          /* worry about transtions too? */
 157                        bCorr;                           /* correlation function ? */
 158     real         t, aa, aver, aver2, aversig, fraction; /* fraction trans dihedrals */
 159     double       tfrac = 0;
 160     char         title[256];
 161     real       **dih = NULL;      /* mega array with all dih. angles at all times*/
 162     char         buf[80];
 163     real        *time, *trans_frac, *aver_angle;
 164     t_filenm     fnm[] = {
 165         { efTRX, "-f", NULL,  ffREAD  },
 166         { efNDX, NULL, "angle",  ffREAD  },
 167         { efXVG, "-od", "angdist",  ffWRITE },
 168         { efXVG, "-ov", "angaver",  ffOPTWR },
 169         { efXVG, "-of", "dihfrac",  ffOPTWR },
 170         { efXVG, "-ot", "dihtrans", ffOPTWR },
 171         { efXVG, "-oh", "trhisto",  ffOPTWR },
 172         { efXVG, "-oc", "dihcorr",  ffOPTWR },
 173         { efTRR, "-or", NULL,       ffOPTWR }
 174     };
 175 #define NFILE asize(fnm)
 176     int          npargs;
 177     t_pargs     *ppa;
 178     output_env_t oenv;
 179
 180     npargs = asize(pa);
 181     ppa    = add_acf_pargs(&npargs, pa);
 182     parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
 183                       NFILE, fnm, npargs, ppa, asize(desc), desc, asize(bugs), bugs,
 184                       &oenv);
 185
 186     mult   = 4;
 187     maxang = 360.0;
 188     bRb    = FALSE;
 189     switch (opt[0][0])
 190     {
 191         case 'a':
 192             mult   = 3;
 193             maxang = 180.0;
 194             break;
 195         case 'd':
 196             break;
 197         case 'i':
 198             break;
 199         case 'r':
 200             bRb = TRUE;
 201             break;
 202     }
 203
 204     if (opt2bSet("-or", NFILE, fnm))
 205     {
 206         if (mult != 4)
 207         {
 208             gmx_fatal(FARGS, "Can not combine angles with trn dump");
 209         }
 210         else
 211         {
 212             please_cite(stdout, "Mu2005a");
 213         }
 214     }
 215
 216     /* Calculate bin size */
 217     maxangstat = (int)(maxang/binwidth+0.5);
 218     binwidth   = maxang/maxangstat;
 219
 220     rd_index(ftp2fn(efNDX, NFILE, fnm), 1, &isize, &index, &grpname);
 221     nangles = isize/mult;
 222     if ((isize % mult) != 0)
 223     {
 224         gmx_fatal(FARGS, "number of index elements not multiple of %d, "
 225                   "these can not be %s\n",
 226                   mult, (mult == 3) ? "angle triplets" : "dihedral quadruplets");
 227     }
 228
 229
 230     /* Check whether specific analysis has to be performed */
 231     bCorr  = opt2bSet("-oc", NFILE, fnm);
 232     bAver  = opt2bSet("-ov", NFILE, fnm);
 233     bTrans = opt2bSet("-ot", NFILE, fnm);
 234     bFrac  = opt2bSet("-of", NFILE, fnm);
 235     if (bTrans && opt[0][0] != 'd')
 236     {
 237         fprintf(stderr, "Option -ot should only accompany -type dihedral. Disabling -ot.\n");
 238         bTrans = FALSE;
 239     }
 240
 241     if (bChandler && !bCorr)
 242     {
 243         bCorr = TRUE;
 244     }
 245
 246     if (bFrac && !bRb)
 247     {
 248         fprintf(stderr, "Warning:"
 249                 " calculating fractions as defined in this program\n"
 250                 "makes sense for Ryckaert Bellemans dihs. only. Ignoring -of\n\n");
 251         bFrac = FALSE;
 252     }
 253
 254     if ( (bTrans || bFrac || bCorr) && mult == 3)
 255     {
 256         gmx_fatal(FARGS, "Can only do transition, fraction or correlation\n"
 257                   "on dihedrals. Select -d\n");
 258     }
 259
 260     /*
 261      * We need to know the nr of frames so we can allocate memory for an array
 262      * with all dihedral angles at all timesteps. Works for me.
 263      */
 264     if (bTrans || bCorr  || bALL || opt2bSet("-or", NFILE, fnm))
 265     {
 266         snew(dih, nangles);
 267     }
 268
 269     snew(angstat, maxangstat);
 270
 271     read_ang_dih(ftp2fn(efTRX, NFILE, fnm), (mult == 3),
 272                  bALL || bCorr || bTrans || opt2bSet("-or", NFILE, fnm),
 273                  bRb, bPBC, maxangstat, angstat,
 274                  &nframes, &time, isize, index, &trans_frac, &aver_angle, dih,
 275                  oenv);
 276
 277     dt = (time[nframes-1]-time[0])/(nframes-1);
 278
 279     if (bAver)
 280     {
 281         sprintf(title, "Average Angle: %s", grpname);
 282         out = xvgropen(opt2fn("-ov", NFILE, fnm),
 283                        title, "Time (ps)", "Angle (degrees)", oenv);
 284         for (i = 0; (i < nframes); i++)
 285         {
 286             fprintf(out, "%10.5f  %8.3f", time[i], aver_angle[i]*RAD2DEG);
 287             if (bALL)
 288             {
 289                 for (j = 0; (j < nangles); j++)
 290                 {
 291                     if (bPBC)
 292                     {
 293                         real dd = dih[j][i];
 294                         fprintf(out, "  %8.3f", atan2(sin(dd), cos(dd))*RAD2DEG);
 295                     }
 296                     else
 297                     {
 298                         fprintf(out, "  %8.3f", dih[j][i]*RAD2DEG);
 299                     }
 300                 }
 301             }
 302             fprintf(out, "\n");
 303         }
 304         ffclose(out);
 305     }
 306     if (opt2bSet("-or", NFILE, fnm))
 307     {
 308         dump_dih_trn(nframes, nangles, dih, opt2fn("-or", NFILE, fnm), time);
 309     }
 310
 311     if (bFrac)
 312     {
 313         sprintf(title, "Trans fraction: %s", grpname);
 314         out = xvgropen(opt2fn("-of", NFILE, fnm),
 315                        title, "Time (ps)", "Fraction", oenv);
 316         tfrac = 0.0;
 317         for (i = 0; (i < nframes); i++)
 318         {
 319             fprintf(out, "%10.5f  %10.3f\n", time[i], trans_frac[i]);
 320             tfrac += trans_frac[i];
 321         }
 322         ffclose(out);
 323
 324         tfrac /= nframes;
 325         fprintf(stderr, "Average trans fraction: %g\n", tfrac);
 326     }
 327     sfree(trans_frac);
 328
 329     if (bTrans)
 330     {
 331         ana_dih_trans(opt2fn("-ot", NFILE, fnm), opt2fn("-oh", NFILE, fnm),
 332                       dih, nframes, nangles, grpname, time, bRb, oenv);
 333     }
 334
 335     if (bCorr)
 336     {
 337         /* Autocorrelation function */
 338         if (nframes < 2)
 339         {
 340             fprintf(stderr, "Not enough frames for correlation function\n");
 341         }
 342         else
 343         {
 344
 345             if (bChandler)
 346             {
 347                 real     dval, sixty = DEG2RAD*60;
 348                 gmx_bool bTest;
 349
 350                 for (i = 0; (i < nangles); i++)
 351                 {
 352                     for (j = 0; (j < nframes); j++)
 353                     {
 354                         dval = dih[i][j];
 355                         if (bRb)
 356                         {
 357                             bTest = (dval > -sixty) && (dval < sixty);
 358                         }
 359                         else
 360                         {
 361                             bTest = (dval < -sixty) || (dval > sixty);
 362                         }
 363                         if (bTest)
 364                         {
 365                             dih[i][j] = dval-tfrac;
 366                         }
 367                         else
 368                         {
 369                             dih[i][j] = -tfrac;
 370                         }
 371                     }
 372                 }
 373             }
 374             if (bChandler)
 375             {
 376                 mode = eacNormal;
 377             }
 378             else
 379             {
 380                 mode = eacCos;
 381             }
 382             do_autocorr(opt2fn("-oc", NFILE, fnm), oenv,
 383                         "Dihedral Autocorrelation Function",
 384                         nframes, nangles, dih, dt, mode, bAverCorr);
 385         }
 386     }
 387
 388
 389     /* Determine the non-zero part of the distribution */
 390     for (first = 0; (first < maxangstat-1) && (angstat[first+1] == 0); first++)
 391     {
 392         ;
 393     }
 394     for (last = maxangstat-1; (last > 0) && (angstat[last-1] == 0); last--)
 395     {
 396         ;
 397     }
 398
 399     aver = aver2 = 0;
 400     for (i = 0; (i < nframes); i++)
 401     {
 402         aver  += RAD2DEG*aver_angle[i];
 403         aver2 += sqr(RAD2DEG*aver_angle[i]);
 404     }
 405     aver   /= (real) nframes;
 406     aver2  /= (real) nframes;
 407     aversig = sqrt(aver2-sqr(aver));
 408     printf("Found points in the range from %d to %d (max %d)\n",
 409            first, last, maxangstat);
 410     printf(" < angle >  = %g\n", aver);
 411     printf("< angle^2 > = %g\n", aver2);
 412     printf("Std. Dev.   = %g\n", aversig);
 413
 414     if (mult == 3)
 415     {
 416         sprintf(title, "Angle Distribution: %s", grpname);
 417     }
 418     else
 419     {
 420         sprintf(title, "Dihedral Distribution: %s", grpname);
 421
 422         calc_distribution_props(maxangstat, angstat, -180.0, 0, NULL, &S2);
 423         fprintf(stderr, "Order parameter S^2 = %g\n", S2);
 424     }
 425
 426     bPeriodic = (mult == 4) && (first == 0) && (last == maxangstat-1);
 427
 428     out = xvgropen(opt2fn("-od", NFILE, fnm), title, "Degrees", "", oenv);
 429     if (output_env_get_print_xvgr_codes(oenv))
 430     {
 431         fprintf(out, "@    subtitle \"average angle: %g\\So\\N\"\n", aver);
 432     }
 433     norm_fac = 1.0/(nangles*nframes*binwidth);
 434     if (bPeriodic)
 435     {
 436         maxstat = 0;
 437         for (i = first; (i <= last); i++)
 438         {
 439             maxstat = max(maxstat, angstat[i]*norm_fac);
 440         }
 441         fprintf(out, "@with g0\n");
 442         fprintf(out, "@    world xmin -180\n");
 443         fprintf(out, "@    world xmax  180\n");
 444         fprintf(out, "@    world ymin 0\n");
 445         fprintf(out, "@    world ymax %g\n", maxstat*1.05);
 446         fprintf(out, "@    xaxis  tick major 60\n");
 447         fprintf(out, "@    xaxis  tick minor 30\n");
 448         fprintf(out, "@    yaxis  tick major 0.005\n");
 449         fprintf(out, "@    yaxis  tick minor 0.0025\n");
 450     }
 451     for (i = first; (i <= last); i++)
 452     {
 453         fprintf(out, "%10g  %10f\n", i*binwidth+180.0-maxang, angstat[i]*norm_fac);
 454     }
 455     if (bPeriodic)
 456     {
 457         /* print first bin again as last one */
 458         fprintf(out, "%10g  %10f\n", 180.0, angstat[0]*norm_fac);
 459     }
 460
 461     ffclose(out);
 462
 463     do_view(oenv, opt2fn("-od", NFILE, fnm), "-nxy");
 464     if (bAver)
 465     {
 466         do_view(oenv, opt2fn("-ov", NFILE, fnm), "-nxy");
 467     }
 468
 469     return 0;
 470 }