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