src/gromacs/gmxana/gmx_angle.c

   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, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 #include "config.h"
  38 #include <math.h>
  39 #include <string.h>
  40
  41 #include "gromacs/math/units.h"
  42 #include "typedefs.h"
  43 #include "gromacs/utility/smalloc.h"
  44 #include "gromacs/utility/futil.h"
  45 #include "gromacs/commandline/pargs.h"
  46 #include "copyrite.h"
  47 #include "gromacs/math/vec.h"
  48 #include "gromacs/topology/index.h"
  49 #include "macros.h"
  50 #include "gromacs/utility/fatalerror.h"
  51 #include "gromacs/fileio/xvgr.h"
  52 #include "viewit.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         gmx_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         gmx_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         if (output_env_get_print_xvgr_codes(oenv))
 445         {
 446             fprintf(out, "@with g0\n");
 447             fprintf(out, "@    world xmin -180\n");
 448             fprintf(out, "@    world xmax  180\n");
 449             fprintf(out, "@    world ymin 0\n");
 450             fprintf(out, "@    world ymax %g\n", maxstat*1.05);
 451             fprintf(out, "@    xaxis  tick major 60\n");
 452             fprintf(out, "@    xaxis  tick minor 30\n");
 453             fprintf(out, "@    yaxis  tick major 0.005\n");
 454             fprintf(out, "@    yaxis  tick minor 0.0025\n");
 455         }
 456     }
 457     for (i = first; (i <= last); i++)
 458     {
 459         fprintf(out, "%10g  %10f\n", i*binwidth+180.0-maxang, angstat[i]*norm_fac);
 460     }
 461     if (bPeriodic)
 462     {
 463         /* print first bin again as last one */
 464         fprintf(out, "%10g  %10f\n", 180.0, angstat[0]*norm_fac);
 465     }
 466
 467     gmx_ffclose(out);
 468
 469     do_view(oenv, opt2fn("-od", NFILE, fnm), "-nxy");
 470     if (bAver)
 471     {
 472         do_view(oenv, opt2fn("-ov", NFILE, fnm), "-nxy");
 473     }
 474
 475     return 0;
 476 }