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 "gromacs/fileio/futil.h"
  46 #include "gromacs/commandline/pargs.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 "gromacs/fileio/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         "[THISMODULE] 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, [THISMODULE] 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 [gmx-covar].[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     if (!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         return 0;
 187     }
 188
 189     mult   = 4;
 190     maxang = 360.0;
 191     bRb    = FALSE;
 192     switch (opt[0][0])
 193     {
 194         case 'a':
 195             mult   = 3;
 196             maxang = 180.0;
 197             break;
 198         case 'd':
 199             break;
 200         case 'i':
 201             break;
 202         case 'r':
 203             bRb = TRUE;
 204             break;
 205     }
 206
 207     if (opt2bSet("-or", NFILE, fnm))
 208     {
 209         if (mult != 4)
 210         {
 211             gmx_fatal(FARGS, "Can not combine angles with trn dump");
 212         }
 213         else
 214         {
 215             please_cite(stdout, "Mu2005a");
 216         }
 217     }
 218
 219     /* Calculate bin size */
 220     maxangstat = (int)(maxang/binwidth+0.5);
 221     binwidth   = maxang/maxangstat;
 222
 223     rd_index(ftp2fn(efNDX, NFILE, fnm), 1, &isize, &index, &grpname);
 224     nangles = isize/mult;
 225     if ((isize % mult) != 0)
 226     {
 227         gmx_fatal(FARGS, "number of index elements not multiple of %d, "
 228                   "these can not be %s\n",
 229                   mult, (mult == 3) ? "angle triplets" : "dihedral quadruplets");
 230     }
 231
 232
 233     /* Check whether specific analysis has to be performed */
 234     bCorr  = opt2bSet("-oc", NFILE, fnm);
 235     bAver  = opt2bSet("-ov", NFILE, fnm);
 236     bTrans = opt2bSet("-ot", NFILE, fnm);
 237     bFrac  = opt2bSet("-of", NFILE, fnm);
 238     if (bTrans && opt[0][0] != 'd')
 239     {
 240         fprintf(stderr, "Option -ot should only accompany -type dihedral. Disabling -ot.\n");
 241         bTrans = FALSE;
 242     }
 243
 244     if (bChandler && !bCorr)
 245     {
 246         bCorr = TRUE;
 247     }
 248
 249     if (bFrac && !bRb)
 250     {
 251         fprintf(stderr, "Warning:"
 252                 " calculating fractions as defined in this program\n"
 253                 "makes sense for Ryckaert Bellemans dihs. only. Ignoring -of\n\n");
 254         bFrac = FALSE;
 255     }
 256
 257     if ( (bTrans || bFrac || bCorr) && mult == 3)
 258     {
 259         gmx_fatal(FARGS, "Can only do transition, fraction or correlation\n"
 260                   "on dihedrals. Select -d\n");
 261     }
 262
 263     /*
 264      * We need to know the nr of frames so we can allocate memory for an array
 265      * with all dihedral angles at all timesteps. Works for me.
 266      */
 267     if (bTrans || bCorr  || bALL || opt2bSet("-or", NFILE, fnm))
 268     {
 269         snew(dih, nangles);
 270     }
 271
 272     snew(angstat, maxangstat);
 273
 274     read_ang_dih(ftp2fn(efTRX, NFILE, fnm), (mult == 3),
 275                  bALL || bCorr || bTrans || opt2bSet("-or", NFILE, fnm),
 276                  bRb, bPBC, maxangstat, angstat,
 277                  &nframes, &time, isize, index, &trans_frac, &aver_angle, dih,
 278                  oenv);
 279
 280     dt = (time[nframes-1]-time[0])/(nframes-1);
 281
 282     if (bAver)
 283     {
 284         sprintf(title, "Average Angle: %s", grpname);
 285         out = xvgropen(opt2fn("-ov", NFILE, fnm),
 286                        title, "Time (ps)", "Angle (degrees)", oenv);
 287         for (i = 0; (i < nframes); i++)
 288         {
 289             fprintf(out, "%10.5f  %8.3f", time[i], aver_angle[i]*RAD2DEG);
 290             if (bALL)
 291             {
 292                 for (j = 0; (j < nangles); j++)
 293                 {
 294                     if (bPBC)
 295                     {
 296                         real dd = dih[j][i];
 297                         fprintf(out, "  %8.3f", atan2(sin(dd), cos(dd))*RAD2DEG);
 298                     }
 299                     else
 300                     {
 301                         fprintf(out, "  %8.3f", dih[j][i]*RAD2DEG);
 302                     }
 303                 }
 304             }
 305             fprintf(out, "\n");
 306         }
 307         ffclose(out);
 308     }
 309     if (opt2bSet("-or", NFILE, fnm))
 310     {
 311         dump_dih_trn(nframes, nangles, dih, opt2fn("-or", NFILE, fnm), time);
 312     }
 313
 314     if (bFrac)
 315     {
 316         sprintf(title, "Trans fraction: %s", grpname);
 317         out = xvgropen(opt2fn("-of", NFILE, fnm),
 318                        title, "Time (ps)", "Fraction", oenv);
 319         tfrac = 0.0;
 320         for (i = 0; (i < nframes); i++)
 321         {
 322             fprintf(out, "%10.5f  %10.3f\n", time[i], trans_frac[i]);
 323             tfrac += trans_frac[i];
 324         }
 325         ffclose(out);
 326
 327         tfrac /= nframes;
 328         fprintf(stderr, "Average trans fraction: %g\n", tfrac);
 329     }
 330     sfree(trans_frac);
 331
 332     if (bTrans)
 333     {
 334         ana_dih_trans(opt2fn("-ot", NFILE, fnm), opt2fn("-oh", NFILE, fnm),
 335                       dih, nframes, nangles, grpname, time, bRb, oenv);
 336     }
 337
 338     if (bCorr)
 339     {
 340         /* Autocorrelation function */
 341         if (nframes < 2)
 342         {
 343             fprintf(stderr, "Not enough frames for correlation function\n");
 344         }
 345         else
 346         {
 347
 348             if (bChandler)
 349             {
 350                 real     dval, sixty = DEG2RAD*60;
 351                 gmx_bool bTest;
 352
 353                 for (i = 0; (i < nangles); i++)
 354                 {
 355                     for (j = 0; (j < nframes); j++)
 356                     {
 357                         dval = dih[i][j];
 358                         if (bRb)
 359                         {
 360                             bTest = (dval > -sixty) && (dval < sixty);
 361                         }
 362                         else
 363                         {
 364                             bTest = (dval < -sixty) || (dval > sixty);
 365                         }
 366                         if (bTest)
 367                         {
 368                             dih[i][j] = dval-tfrac;
 369                         }
 370                         else
 371                         {
 372                             dih[i][j] = -tfrac;
 373                         }
 374                     }
 375                 }
 376             }
 377             if (bChandler)
 378             {
 379                 mode = eacNormal;
 380             }
 381             else
 382             {
 383                 mode = eacCos;
 384             }
 385             do_autocorr(opt2fn("-oc", NFILE, fnm), oenv,
 386                         "Dihedral Autocorrelation Function",
 387                         nframes, nangles, dih, dt, mode, bAverCorr);
 388         }
 389     }
 390
 391
 392     /* Determine the non-zero part of the distribution */
 393     for (first = 0; (first < maxangstat-1) && (angstat[first+1] == 0); first++)
 394     {
 395         ;
 396     }
 397     for (last = maxangstat-1; (last > 0) && (angstat[last-1] == 0); last--)
 398     {
 399         ;
 400     }
 401
 402     aver = aver2 = 0;
 403     for (i = 0; (i < nframes); i++)
 404     {
 405         aver  += RAD2DEG*aver_angle[i];
 406         aver2 += sqr(RAD2DEG*aver_angle[i]);
 407     }
 408     aver   /= (real) nframes;
 409     aver2  /= (real) nframes;
 410     aversig = sqrt(aver2-sqr(aver));
 411     printf("Found points in the range from %d to %d (max %d)\n",
 412            first, last, maxangstat);
 413     printf(" < angle >  = %g\n", aver);
 414     printf("< angle^2 > = %g\n", aver2);
 415     printf("Std. Dev.   = %g\n", aversig);
 416
 417     if (mult == 3)
 418     {
 419         sprintf(title, "Angle Distribution: %s", grpname);
 420     }
 421     else
 422     {
 423         sprintf(title, "Dihedral Distribution: %s", grpname);
 424
 425         calc_distribution_props(maxangstat, angstat, -180.0, 0, NULL, &S2);
 426         fprintf(stderr, "Order parameter S^2 = %g\n", S2);
 427     }
 428
 429     bPeriodic = (mult == 4) && (first == 0) && (last == maxangstat-1);
 430
 431     out = xvgropen(opt2fn("-od", NFILE, fnm), title, "Degrees", "", oenv);
 432     if (output_env_get_print_xvgr_codes(oenv))
 433     {
 434         fprintf(out, "@    subtitle \"average angle: %g\\So\\N\"\n", aver);
 435     }
 436     norm_fac = 1.0/(nangles*nframes*binwidth);
 437     if (bPeriodic)
 438     {
 439         maxstat = 0;
 440         for (i = first; (i <= last); i++)
 441         {
 442             maxstat = max(maxstat, angstat[i]*norm_fac);
 443         }
 444         fprintf(out, "@with g0\n");
 445         fprintf(out, "@    world xmin -180\n");
 446         fprintf(out, "@    world xmax  180\n");
 447         fprintf(out, "@    world ymin 0\n");
 448         fprintf(out, "@    world ymax %g\n", maxstat*1.05);
 449         fprintf(out, "@    xaxis  tick major 60\n");
 450         fprintf(out, "@    xaxis  tick minor 30\n");
 451         fprintf(out, "@    yaxis  tick major 0.005\n");
 452         fprintf(out, "@    yaxis  tick minor 0.0025\n");
 453     }
 454     for (i = first; (i <= last); i++)
 455     {
 456         fprintf(out, "%10g  %10f\n", i*binwidth+180.0-maxang, angstat[i]*norm_fac);
 457     }
 458     if (bPeriodic)
 459     {
 460         /* print first bin again as last one */
 461         fprintf(out, "%10g  %10f\n", 180.0, angstat[0]*norm_fac);
 462     }
 463
 464     ffclose(out);
 465
 466     do_view(oenv, opt2fn("-od", NFILE, fnm), "-nxy");
 467     if (bAver)
 468     {
 469         do_view(oenv, opt2fn("-ov", NFILE, fnm), "-nxy");
 470     }
 471
 472     return 0;
 473 }