src/gromacs/gmxana/gmx_angle.cpp

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