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